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The  task  of  quantitatively  describing  physical  properties  of 
solutions,  whether  for  use  in  engineering  design  and  optimization,  or 
for  understanding  natural  processes,  grows  increasingly  difficult  as 
more  accuracy  is  demanded  and  more  complex  systems  are  studied.   To 
meet  these  demands,  engineering  thermodynamics  now  looks  more  thoroughly 
into  the  intermolecular  forces  in  a  solution,  and  rely  more  on  the  the- 
oretical,  rather  than  empirical,  bases  for  their  quantitative  descrip- 
tions.  This  molecular  basis  is  particularly  important  when  "reactive" 
systems — those  systems  where  component  molecules  interact  so  as  to  form 
different  chemical  species — are  considered.   This  work  develops  a  sta- 
tistical mechanical  solution  theory  for  such  systems,  examining  the 
constraints  caused  by  their  "reactive"  nature,  and  a  model  which  may  be 
applied  to  correlate  and  predict  their  thermodynamic  properties. 


xii 


The  thermodynamic  consistency  requirements  of  models  for  "reactive" 
solutions  are  expressed.   These  are  used  to  formulate  a  procedure,  using 
matrix  projections,  which  builds  these  requirements  into  the  solution 
model.   The  statistical  mechanics  of  such  a  solution  in  the  grand  canon- 
ical ensemble  is  explored  in  light  of  the  "reactive"  constraints  and  the 
consistency  requirements.   By  using  the  projection  procedure  mentioned 
above,  the  derivative,  or  fluctuation,  properties  of  the  components  in 
a  general  "reactive"  system  are  related  to  equilibrium  relations  and 
molecular  correlation  functions. 

The  case  of  total  dissociations  as  the  only  "reactions"  in  the 
solution  is  explored,  with  specific  applications  to  strong  electrolytes. 
First,  the  singularities  of  the  long  range  interactions  are  shown  to  be 
removed  by  the  projection  procedure.   Then  a  model  is  formulated  based 
on  a  composite  of  hard-sphere  and  charge  effects.   This  model,  though 
presently  not  complete,  shows  considerable  promise  in  correlating  the 
properties  of  salt  solutions.   An  experimental  apparatus  for  volumetric 
behavior  of  liquids  was  completed  and  methods  are  described  for  using 
data  from  it  to  evaluate  the  parameters  in  the  solution  model,  as  well 
as  develop  a  better  understanding  of  density  effects  in  electrolyte 
solutions. 
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CHAPTER  1 
INTRODUCTION 


Design  and  optimization  of  industrial  processes  and  equipment, 
as  well  as  quantitative  understanding  of  natural  or  biological  pro- 
cesses, require  quantitative  descriptions  of  the  physical  properties 
of  the  components  of  those  systems.   A  very  important  need  is  the  cor- 
relation and  prediction  of  P-V-T-x.  relationships  and  phase  equilibria 
in  multicomponent  liquid  systems.   Meeting  this  need  can  be  quite  dif- 
ficult, especially  if  there  are  complicated  intermolecular  interactions 
or  tight  constraints  on  accuracy.   The  situation  is  complicated  further 
in  "reactive"  systems,  where  the  relative  amounts  of  the  species  in  a 
given  system  are  constrained  by  equilibrium  conditions  or  stoichiometry . 
These  "reactive"  systems  include,  in  general,  solvation,  dissociation 
and  association  —  including  processes  such  as  micellization  or  ioniza- 
tion of  salts.   A  large  number  of  real  systems  exhibit  this  "reactive" 
nature,  although  it  may  proceed  to  such  a  small  degree  that  it  is 
negligeable.  In  other  situations,  however,  the  additional  species  and 
constraint  relations  contribute  appreciably  to  the  solution  properties. 
Although  specific  types  of  systems  have  been  treated  previously  (e.g. 
Renon  and  Prausnitz,  1967),  there  has  not  been  a  general  approach  to 
these  systems  which  could  serve  as  a  guideline  for  the  consistent  use 
of  solution  models. 


Solutions  involving  electrolytes  are  a  class  of  "reactive"  systems 
which  are  of  particular  interest.   Electrolytes  appear  in  geological 
liquid  systems  (oil  reservoirs,  geothermal  energy  sources),  in  systems 
of  biological  interest  (e.g.  in  active  transport  in  membranes  and  in 
using  salts  to  separate  proteins),  and  have  important  applications  to 
industrial  systems  (extraction,  distillation  and  reverse  osmosis). 
These  systems  have  an  additional  complication,  though,  because  of  the 
presence  of  long-range  coulombic  interactions  in  addition  to  equili- 
brium and  stoichiometric  constraints. 

Statistical  mechanical  approaches  have  been  successful  for  many 
non-"reactive"  systems,  as  well  as  some  specific  "reactive"  systems 
(Friedman  and  Ramanathan,  1970) .   A  recent  approach  using  fluctuation 
solution  theory  (Mathias  and  O'Connell,  1980a, b)  has  shown  application 
for  gases  dissolved  in  various  solvents  at  high  pressures.   This  pro- 
cedure is  based  on  molecular  correlation  functions  and  works  particu- 
larly well  in  systems  which  contain  components  of  very  different  na- 
tures.  Hence,  it  should  be  useful  for  salt  solutions. 

This  work  aims,  firstly,  to  treat  "reactive"  systems  in  a  manner 
that  is  applicable  to  any  system,  and  so  to  present  a  general  and 
unifying  thermodynamic  framework  in  which  correlations  can  be  consis- 
tently developed.   This  framework  is  then  applied  to  the  above  statis- 
tical mechanical  method  to  develop  a  formalism  for  expressing  the 
properties  of  "reactive"  solutions  in  terms  of  these  molecular  correla- 
tion functions.   Finally,  this  formalism  is  applied  to  solutions  of 
strong  electrolytes  and  a  solution  model  is  examined. 


Chapter  two  develops  the  classical  thermodynamics  of  "reactive" 
systems.   There,  the  "reactive"  nature  is  defined  mathematically  in 
terms  of  the  material  balance  (stoichiometric)  relations  and  intensive 
extents  of  reactions.   The  total  amounts  of  chemical  species  in  a  sys- 
tem are  related  to  the  amounts  of  components  (where  the  components  each 
represent  a  degree  of  freedom  in  composition)  by  a  matrix  relation  which 
is,  in  effect,  a  projection  of  the  species  space  into  the  component  sub- 
space.   The  projection  operators  are  completely  defined  by  stoichiometry 
and  equilibria,  and  are  used  to  relate  solution  properties  in  terms  of 
species   to  the  required  solution  properties  in  terms  of  components. 
These  relations  also  show  restrictions  which  must  be  placed  on  solution 
models  for  such  systems,  and  two  specific  examples  are  detailed. 

Chapter  three  begins  with  the  definition  of  a  grand  ensemble  with 
overall  constraints  on  the  mole  numbers  of  the  species.   It  then  derives 
a  solution  theory  based  on  total  and  direct  correlation  functions, 
between  species  that  yields  component  properties.   This  functionality 
is  important  because,  although  the  ensemble-averaged  species  composition 
vector  must  lie  in  the  component  subspace,  this  is  not  necessary  for 
the  species  composition  vectors  for  any  of  the  systems  in  the  ensemble. 
The  chapter  then  gives  the  component  density  derivatives  of  the  com- 
ponent activities  in  terms  of  the  reaction  extents  and  the  integrals 
of  the  species  direct  correlation  function  integrals.   Lastly,  a  de- 
tailed analysis  of  the  matrix  equations  for  a  system  which  undergoes 
complete  dissociations  is  given,  and  the  matrices  are  shown  to  col- 
lapse to  a  very  simple  form. 


Chapter  four  considers  a  specific  case  of  complete  dissociation — 
a  system  of  strong  electrolytes  dissolved  in  a  non-  "reactive"  solvent. 
First,  it  reformulates  the  results  of  chapter  three  using  dimensionally 
smaller  matrices,  explicitly  accounting  for  the  absence  of  the  undis- 
sociated  salts.   In  this  formulation,  it  is  shown  that  the  long-range 
portions  of  the  direct  correlation  function  integrals  (which  diverge 
in  integration)  are  exactly  removed  by  the  projections.   Isothermal 
compressibilities,  partial  molar  volumes  and  composition  derivatives 
of  the  activity  coefficient  are  given  in  terms  of  special  integrals 
of  the  ion-ion,  ion-solvent,  and  solvent-solvent  direct  correlation 
functions,  and  the  Debye-Hueckel  limits  of  these  properties  are  given. 
Two  specific  cases  are  examined,  the  simplest  case  of  one  salt  and 
one  solvent  and  the  case  of  two  salts  with  a  common  anion.   The  latter 
example  gives  the  Harned  coefficients  in  terms  of  the  correlation  func- 
tion integrals. 

Chapter  five  deals  with  the  development  of  a  model  for  the  elec- 
trolyte solutions.   Results  of  previous  theory  and  experimental  results 
are  used  to  determine  the  type  of  model  to  be  used.   The  direct  cor- 
relation function  between  molecules  i  and  j  is  given: 


C.j(r)  =  C.^^(r)  +   {C..(r)  -  C.'^^(r)} 


where  the  terms  superscripted  hs  are  the  direct  correlation  functions 
for  a  system  of  hard  spheres.   The  bracketed  term,  which  corresponds 
primarily  to  electrostatic  effects,  is  approximated  by  a  composite 


expression,  based  on  a  low-density  expansion  such  as  the  mean  spherical 
model,  with  correction  to  the  proper  infinite  dilution  properties  in- 
cluding hydration  of  ions.   The  model  is  formulated  in  terms  of  several 
universal  (not  salt  dependent)  constants,  and  functions  of  ionic  size, 
charge   and  density.   Some  rough,  preliminary  calculations  are  given. 

Chapter  six  examines  an  experimental  method  by  which  volumetric 
behavior  of  salt  solutions  can  be  determined  over  a  sizeable  range  of 
of  pressures  and  temperatures.   The  range  of  operation  allows  the  ef- 
fects of  total  density  and  composition  on  free  energies  to  be  examined 
individually  by  allowing  available  activity  data  at  low  pressures  to 
be  extended  to  higher  pressures   (at  constant  density) .   It  also  pro- 
vides a  basis  for  determining  the  model  parameters  for  the  electrolyte 
direct  correlation  functions.   The  design  and  operation  of  this  ex- 
periment are  described. 

Finally,  further  applications  of  the  results  of  this  work,  as 
well  as  suggestions  for   extensions,  are  given. 


CHAPTER  2 
SOLUTION  THERMODYNAI-IICS  FOR  "REACTIVE"  COMPONENTS 


Introduction 

Fundamental  to  the  development  of  thermodynamic  properties  is  the 
mathematical  rigor  that  yields  basic  relationships  such  as  Maxwell  rela- 
tions (Van  Ness,  1964),  identification  of  partial  derivatives  (e.g., 
heat  capacities),  and  Legendre  Transformations  (Beegle  et  al . ,  197A) . 
Provided  the  requirements  of  the  specification  of  thermodynamic  state 
are  met,  natural  phenomena  appear  to  follow  all  the  mathematical  conse- 
quences of  the  First  and  Second  Laws  and  subsequent  definitions,  even 
for  systems  of  many  components  and  modes  of  interaction  with  the  sur- 
roundings (Redlich,  1968). 

The  most  complex  systems  to  deal  with  are  those  involving  chemical 
reactions  because  the  compositions  of  the  species  are  connected  at  equi- 
librium. In  regard  to  this  fact,  the  purpose  of  this  chapter  is  to 
display  a  mathematical  formalism  for  expressing  the  solution  thermo- 
dynamic properties  of  mixtures  made  of  components  which  can  react  to 
form  additional  species.   Their  concentrations  can  be  related  through 
equilibrium  constants  or  there  can  be   complete  dissociation.   The 
particular  situations  of  current  interest  include  electrolytes,  both 
strong  and  weak,  the  "chemical  theory"  of  solutions  (Prausnitz,  1969) 
and  the  "solution  of  groups"  method  for  activity  coefficients  (Derr  and 
Deal,  1968,  1969,  1973;  Fredenslund  et  al.,  1977a).  Our  motivation  is  to' 

6 


provide  a  degree  of  unification  not  heretofore  attempted  so  that  deri- 
vations can  be  streamlined  and  any  potential  inconsistencies  can  be 
avoided,  particularly  in  the  development  of  molecular  and  statistical 
thermodynamic  theories  and  correlations.   The  results  obtained  for  the 
first  two  cases  are  the  same  as  found  previously.   They  have  been  ex- 
tended to  a  more  complete  set  of  thermodynamic  properties.   In  the 
solution  of  groups  method  we  have  found  precisely  how  analytical  solu- 
tion of  groups  methods  such  as  UNIFAC  (Fredenslund  et  al.,  1977b)  yield 
different  results  from  the  molecular  equation  upon  which  they  are  based 
(e.g.,  UNIFAC  from  UNIQUAC,  Abrams  and  Prausnitz,  1975). 

The  results  of  this  analysis  form  the  framework  for  the  more 
detailed  analysis  of  the  statistical  mechanics  of  solutions  of  "reac- 
tive" components,  and  in  particular  for  the  modelling  of  properties  of 
solutions  of  strong  electrolytes,  which  are  the  subjects  of  subsequent 
chapters. 

Mathematical  Basis 
The  mathematical  situation  is  one  in  which  the  thermodynamic 
functions  of  interest  (enthalpy,  entropy,  activities,  etc.)  are  most 
conveniently  expressed  in  terms  of  a  set  of  variables  related  to  the 
number  of  moles  of  species,  N,  or  their  mole  fractions,  x,    but  actu- 
ally must  be  functions  of  the  set  of  independent  values  related  to 

the  number  of  moles  of  the  original  components,  N  ,  or  their  mole 

o 

fractions,  x  .   The  fact  that  the  number  of  independent  variables  in 
o 

the  former  case  is  different  from  the  number  in  the  latter  case  re- 
quires that  the  functions  be  projected  from  one  space  to  another. 


Figure  1  Is  a  graphic  representation  of  the  projection  process. 
Here,  we  treat  the  case  of  three  species,  with  two  independent  mole 
fractions,  in  a  two  component  system  with  only  one  independent  mole 
fraction.   The  surface  represents  the  values  of  the  thermodynamic 
property,  M,  for  all  values  of  the  species  fractions  satisfying 


^  X.  =  1 
.   1 

1 


However,  the  "allowable"  compositions  of  species  are  controlled  by 
the  stoichiometry  of  the  "reaction"  system.   These  values  are  con- 
strained to  the  projection  of  the  species  fractions,  which  here  is  a 
line  in  the  base  plane.   The  attainable  values  of  the  system  thermo- 
dynamic property,  M^,  are  the  projection  of  M,  the  curve  formed  by 
the  intersection  of  its  surface  with  the  vertical  plane  containing  the 
line  of  "allowed"  species  fractions.   In  multicomponent-multispecies 
systems  the  system  thermodynamic  property  values  will  be  found  in  the 
intersection  of  the  property  surface  in  species  space  and  a  hypersur- 
face  for  the  allowable  species  fractions.   Our  purpose  is  to  establish 
a  general  procedure  to  obtain  this  intersection.   The  mathematical  op- 
erators for  this  procedure  in  the  present  case  can  be  straightforwardly 
defined  and  are  easily  used  as  standard  matrix  manipulations."^ 

This  method  has  been  implicitly  used  in  a  wide  variety  of  sit- 
uations such  as  electrolyte  solutes,  the  "chemical  theory"  of  solu- 
tions and  "analytical  solution  of  groups"  correlations.   While  the 


1 
In  certain  statistical  thermodynamic  developments  the  method  is 
more  complex  but  represents  the  only  technique  we  know  to  remove  cer- 
tain apparent  divergences  (Perry  et  al.  ,  19vS0a) 


Figure,!.   Projection  of  a  thGrmodynnmic  property  surfaqe. 
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number  of  "new"  results  is  not  large,  our  objective  is  to  unify  appar- 
ently diverse  approaches  and  present  a  single  procedure  for  use  in  all 
problems  where  the  properties  of  a  solution  are  to  be  obtained  from 
model  expressions  for  the  properties  of  the  species. 

Construction  of  Projectors —  Material  Balance 

We  will  begin  by  examining  a  closed  system  at  fixed  T  and  P  made 

from  the  mixing  of  n   components,  with  the  mole  numbers  of  the  com- 

o 

ponents  given  by  the  vector  N  .   This  system  may  now  undergo  successive 
associations,  dissociations  or  reactions  leading  to  an  equilibrium  sys- 
tem of  n  species,  with  mole  numbers  given  by  the  vector  N.   In  order  to 
relate  the  species  vector  N  to  the  component  vector  N  we  will  use  a 
process  which  is  chosen  for  mathematical  convenience  and  is  not  intended 
to  reflect  the  chemical  evolution  of  the  system.   The  general  develop- 
ment is  given  in  Appendix  A;  two  complete  examples  are  given  in  Appendix 
B. 

The  general  result  is 

N  =  WN  (2-1) 

o 

where  W  (defined  in  Appendix  A)  is  a  matrix  involving  the  stoichiomet- 
ric coefficients  and  extents  of  reaction  of  the  independent  reactions 
which  change  components  into  species. 

The  importance  of  this  relationship  lies  in  the  fact  that  N   has 

~o 

n   degrees  of  freedom,  and  since  the  extents  of  reaction  are  functions 
of  N   only  (at  constant  T  and  P) ,  N  also  has  only  n   degrees  of  free- 
dom.  However,  it  is  in  a  space  with  n  degrees  of  freedom.   Thus  W  is 
a  projector  from  a  space  of  dimension  n   into  a  space  of  dimension  n. 
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As  shown  below,  use  of  this  and  related  matrices  allows  us  to  project, 
without  any  loss  of  information,  expressions  for  thermodynamic  proper- 
ties into  the  desired   n^-  space  that  are  naturally  written  in  n-space. 
It  must  be  emphasized  that  the  apparently  simple  relationship  of  equa- 
tion (2-1)  contains  a  considerable  number  of  subleties.   The  next  few 
sections  show  the  ramifications  that  this  relation  has  on  the  thermody- 
namics of  solutions.   An  example  involving  2  components,  4  species,  and 
2  reactions  is  given  in  Appendix  B. 

Properties  of  Solutions 
Partial  Molar  Properties 

We  will  now  see  how  solution  properties,  especially  partial  molar 
properties,  are  projected  using  the  W  matrix  defined  above.   First, 
assuming  our  mixture  to  have  reached  equilibrium  and  for  T  and  P  to  be 
constant,  we  know 


dG)     =  0  where   G  =  n'^P  (2-2) 

T,P 

3G 


and  y .  =  yj^ 


1 


T,P,N 


Since  G  is  really  only  a  function  of  the  N  with  T  and  P  constant, 
we  can  use  Euler's  theorem  to  also  write 


T  ^C 

G  =  N  y      with      u    ="5^7 

— cr^o  ^oa    3N 


oa- 


(2-3) 

'    oBy^a 
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Reaction  equilibrium  can  be  expressed  as 


T 
V  y  =  0 


(2-4) 


where  v  is  an  n  x  r  matrix  with  (v).,  being  the  stoichiometric  coef- 

—  1 K 

ficient  of  species  i  in  reaction  k  and  r  being  the  total  number  of 
independent  reactions  (see  Appendix  H) .   Then,  after  some  manipulation 
(see  Appendix  A) ,  we  obtain 


T 
W  y  =  y 
^0 


(2-5) 


A  consequence  of  this  is  the  well-known  result  that  for  any  component 
a  which  is  also  found  as  a  species  in  the  solution, 


y    =  y 

a         oa 


(2-6) 


or 


9G 
9N 


a 


3G 

3N 


T,P,N 


oa 


j/ot 


T,P,N 


O&^OL 


(2-7) 


Then,  since  all  the  partial  molar  quantities  of  interest  are  linear 

combinations  of  y   (or  y)  and  its  derivatives,  all  of  the  other  partial 
~-o     ~ 

molar  properties  of  components  are  also  unchanged  by  the  reactions. 
That  is,  for  a  general  partial  molar  quantity  (M  =  volume,  energy, 
entropy,  log  of  fugacity,  etc.)  we  have: 


M   =  M 
oa    a 


a  =  1, . .  .  ,n 


(2-8) 
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Equation  (2-7)  is  surprisingly  simple  since  the  derivative  in  n  -space 

o   ^ 

is  complex.   However,  as  long  as  M  is  a  linear  combination  of  G  and  its 
derivatives,  equation  (2-8)  holds  because  of  the  linearity  of  (2-2), 
This  is  valuable  in  the  chemical  theory  of  solutions. 

Further,  these  same  linear  relations,  when  applied  to  equation 
(2-4)  also  yield 


v'^M  =  0  (2-9) 


Mixing  Properties 

For  mixing  properties,  the  differences  in  reference  states  of  the 
components  and  species  appear.  Looking,  for  example,  at  the  change  in 
Gibbs  free  energy  on  mixing, 


AG' 


MIX  =  N^AG^IX  =  N^y  -  NTy°  =   RT  N^  ilnt  (2-10) 


while 


AG^^^  =  N^  AG^^^  =  nV   -  N^°  -=   RT  N^  £na  (2-11) 

o     ~o  —      ^o^o    — o— o      — o  — — o 


where  the  relationship  between  N  \i   and  N  \i        depends  upon  the  choice 

o  o 

of  component  reference  state  composition.   From  equation  (2-5),  equa- 
tions (2-10)  and  (2-11)  and  the  Gibbs-Duhem  equation 


T7;MIX   ,T— tMIX    ,  o   ,T  o,  .„  ^^. 

AG    =  W  AG    -  (y  -  W  VI  )  (2-12a) 

— o "o 
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In  terms  of  the  initial  components,  the  two  forms  for  the  chemical  poten- 
tial are 


(u  )  E  ij   =  (y°)   +  R  T  {(£nx  )   +  (iny   )  }  (2-13a) 

— o  a    oa    ~o  a  o  a      o  a 


'"^u")   +  R  T  {(W^£nx)   +  (W^£nT)  }  (2-13b) 


(W  y  )   +  R  T  {(W  £nx)   +  (W  Uny)] 
—  —    a  a ot 


Equation  (2-13a)  is  the  normal  form  for  y^^.   Equation  (2-13b)  is  the 
general  expression  for  this  quantity  in  terms  of  the  species  variables: 
it  is  applicable  to  all  cases  regardless  of  the  reference  states  chosen. 

The  same  relationships  exist  for  all  partial  molar  properties 
since  they  can  all  be  found  by  linear  operations  (including  differen- 
tiation and  multiplicative  functions)  on  the  chemical  potentials.   Thus 
for 


M.  =  f(y.) 
1      1 


then 


M^ix  =  w^aS  ^«x  -  f  (y°  -  wV°)  (2-12b) 


-Q  O 

Figure  2  shows  a  representation  of  the  projection  method  for 

ah'    when  the  reference  state  composition  for  component  1  is  pure, 

X   =1,  and  "infinite  dilution"  (or  a  hypothetical  state  of  pure  com- 
ol 

ponent  whose  fugacity  is  Henry's  constant)  for  the  others.   Such  a 
choice  is  made  with  electrolytes  and  supercritical  substances.   Here, 
it  has  been  assumed  that  species  2  and  3  are  contained  only  in  component 


15 
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Fip;ure.2.   Projection  v/itli  "infinite  dilution"  reference  states, 
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2,  so  an  ideal  solution  occurs  when  x  „  =  x„  =  x^  =  0.   The  value  of 
'  oz    2     J 

—  — oo 

M   in  the  reference  state  is  M   which  corresponds  to  no  particular 
o2  o2 

pure  species  value.   It  is  merely  the  extrapolation  to  (hypothetical) 

pure  component  2  (x°„  =1)  of  the  tangent  to  the  line  of  the  projection 

at  X  , .   The  line  of  extrapolation  is  the  ideal  solution  line  which 
ol 

becomes  horizontal  in  the  plot  of  Am    vs.  x  ^.   Note  that  we  have 

o        o2 

terminated  the  plane  at  an  arbitrary  point  where  the  physical  state 
might  change  due  to  a  solid  phase  or  a  critical  point  being  encoun- 
tered.  The  last  term  in  equation  (2-12b)  would  be  the  difference  be- 
tween the  line  shown  and  the  line  of  intersection  of  the  plane  connec- 
ting all  of  the  M.  =  lim  M  and  the  (vertical)  stoichiometric  plane. 
^   x.-^l 

X 

These  two  lines  need  not  be  parallel. 

On  the  other  hand,  Figure  3  indicates  the  results  for  a  case  where 
the  reference  states  for  the  components  are  pure  component.   Thus 
x°   =  1  and  x°2  -    1-   This  graph  of  AM^    vs.  x  2  is  more  likely  to 
show  a  sign  change  even  though  the  surface  in  species  space  might  not 
indicate  such  behavior  when  examined  by  itself.   Slight  undulations  of 
the  surface  can  yield  sigmoidal  projections. 

Among  the  many  ways  of  dealing  with  equations  (2-13),  several 
common  ones  have  been  chosen.   The  method  of  Prigogine  and  Defay  (1954) 
has  the  same  reference  state  for  any  initial  component  which  may  also 
be  present  as  a  species  in  the  solution.   That  is 


U   =  y  a  =  l,...,n  (2-14a) 
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Figure  3.   Projection  with  pure  component  reference  states. 
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o    To 
y  =  L  y  (2-14b) 

~o 


This  is  most  useful  for  solutions  where  the  components  are  considered 
pure,  but  react  when  mixed  with  other  components  (e.g.  solvation). 
It  is  this  case,  and  only  this  one,  for  which  the  most  common  relation 
holds 


Inx     Y   =  ^nx  Y  (2-14c) 

oa  oa     a 'a 


T 
A  second  choice  is  the  case  where  W   is  independent  of  composi- 
tion.  It  is  normally  used  when  complete  reactions  are  assumed.   Ex- 
amples include  strong  electrolytes  and  solutions  of  groups,  which  are 
discussed  in  detail  later.   Here,  it  is  assumed  that 


o    To 
y   =  W  y  (2-15a) 

o 


with 


T       T 
Inx     +  iny     =   W  £nx  +  W  Iny  (2-15b) 

o   ^ 


In  such  a  case  W  is  a  constant  and  there  can  be  no  successive  reactions 
so  that  the  upper  n  xn   elements  of  W  are  either  unity  (unreacting 


solvents)  or  either  zero  or  ratios  of  stoichiometric  coefficients  ( 


com- 


pletely reacting  components)  and  the  lower  (n-n  )xn   elements  form  the 
matrix  v.   It  is  thus  more  convenient  to  reformulate  the  expressions  in 
terms  of  a  matrix  which  has  only  nonzero  rows  as  discussed  below. 
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Finally,  there  are  cases  in  which  the  species  standard  states  bear 
no  relation  to  the  component  standard  states.  An  example  of  this  is  an 
associating  substance  which  even  when  "pure"  is  made  up  of  a  number  of 
species  (Figure  3  is  an  example  of  such  a  system) .  Then  the  full  equa- 
tions (2-13)  must  be  used.  The  technique  here  is  to  subtract  equation 
(2-13b)  for  the  component  reference  state  from  the  expression  at  the 
solution  conditions.   For  equation  (2-13b)  this  is  written  as 


M     -   M°   =   l/{y°  +  AG^'''}  -  lim  wT{y°  +  AG^^^} 
-o   -o   -  -    —       ^o  -  -    — 

-(a) 

T  -  MIX        T  -^ix  T        To 

W  AG     -  lim  W  AG^    +  [  W  -  lim  W  ]  y °  (2-16a) 

-°a)  -°(a) 


where  x,  .  is  the  vector  of  species  mole  fractions  at  the  reference 
-(a) 

state  of  component  a.   In  the  case  of  systems  where  the  components  are 

— j^XX 
present  as  species  and  the  model  for  AG    includes  all  the  contribu- 

A^MIX    ,  .   , 
tions  to  AG     ,  this  becomes 
o 


y  -  y°  =  Jag^IX  _  lT  lin,  ag^IX  (2-16b) 

-(a) 


-o   -o   -—     -  ^o 


If  the  standard  state  for  species  is  chosen  such  that  (x  ) .  =1  and 

—  1 

the  phase  is  the  same  as  that  of  the  solution,  the  last  term  in  equa- 
tion (2-14)  removes  any  composition  independent  terms  arising  from 
the  differences  in  phase  or  structure  in  the  expression  chosen  for 
AG     .   An  example  of  this  procedure  is  given  in  Appendix  B. 
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Often,  the  reactions  are  described  in  terms  of  equilibrium  con- 
stants which  are  defined  as 


JlnK  =  -v'^y°/RT  (2-17a) 


which,  with  equations  (2-4)  and  (2-10) ,  yields 


£nK  =  v^^lna  (2-17b) 


Equations  (2-1)  and  (2-17a)  allow  us  to  rearrange  (2-12)  to  a 
relation  explicit  in  the  extents  of  reaction  and  equilibrium  constants. 
The  general  result  is  given  in  Appendix  A,   Equation  (A-11)  for  one  set 
of  reactions  (no  sequential  reactions)  reduces  to 


W  Ag"'   =  AG^'^^  +  Rir^nK  (2-18a) 

o 


This  relation,  or  equation  (A-11),  can  be  used  in  place  of  equation 
(2-12a)  in  any  of  the  above  derivations,  including  an  alternative  to 
equation  (2-16a) . 

The  other  mixing  properties  have  similar  relationships.   For  one 
set  of  reactions,  we  can  in  general  write 


M.  =  f(y.) 
1      1 
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where  f  is  a  linear  operator  (including  derivatives  and  multiplicative 
functions).   Then 


W^AM^^^  =  AM^^^  +  R^   L(T&nK)  (2-18b) 


For  example,  using  standard  operations  on  equation  (2-15)  for 
definitions  of  y.  not  functions  of  pressure 

W^AH^^^  =  AH^^^  +  m^   ^  (2-18C) 
— o      -  dl/T 

^T^IX  ^  ^IX  (2-18d) 

o 

T_^IX   —MIX     T  d(TAnK) 
^i^IX  ^  ^gfUA  _  j^^i  (2-18e) 

— o      -    dT 


and  so  on.   Thus  equation  (2-18c)  gives  the  familiar  relation  for  the 
heat  effect  associated  with  mixing  the  components  expressed  in  terms 
of  the  activities  of  the  species  and  the  equilibrium  constant 


MIX    T ^MIX     T  T  ^ii/T 

Q  =  AH    =  N  Ah"^""  =  N  W 

O       ~0   O         O     di/r 


d£nK 

-  RN  -—7- 

P,N    -°  ^1/T 


T 
RN 


ajlny 

-R  n'  C   -t-t-  (2-19) 

91/T-  P,N   ~°    ^^/"^ 


^  dJlnK 
,T  >-T 


If  the  temperature  variation  of  the  activity  coefficients  is  ignored, 
the  result  is  that  usually  found  for  reactions  in  the  standard  state. 
This  relation  is  also  appropriate  for  the  chemical  theory  of  solutions, 
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Excess  Properties 

Excess  properties  differ  from  mixing  properties  because  an  addition- 

F   F   F 
al  composition  term  appears  in  some  excess  properties  (G  ,A  ,S  ).   Using 

the  usual  notation  we  have  for  component  values: 


G^/RT  =  N'^g.na   -  ^Inx    =     N  Zny 
°       o   o     o    o    o   o 


P      ^   3T£na 

S^/R  =  /      ° 
o     -o   3T 


P,N 


+N  £nx    etc. 
~o   o 


For  the  species  values,  however,  we  have 


G^/RT  =  N  £na  -  N^^^nxE  U^Zny 


E,     T 
S  /R  =  N 


9T£n£ 


-   8T 


+  N'^£,nx 


P,N 
— o 


Use  of  equation  (2-5)  on  the  above  relations  gives 


G^  =  ]/g^  -  Klilnx     -   W^^nx)  -  (y°  -  W^y°) 
-o o '-^o  -  - 


(2-20a) 


=  W^AG^^^  -  RT£nx  -  (m°  -  W^y°) 
o    — o ■ 


(2-20b) 


For  the  case  of  equations  (2-14) ,  the  activity  coefficients  become 


iny       =   £nY  +  S,n(x  /x  ) 
oa     a      a  oa 


(2-21a) 
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For  the  case  of  equations  (2-15)  the  relation  is 


Iny       =  (W^Jlny)   +  (w'^ilnx)   -  Inx  (2-21b) 

oa a a     oa 


For  the  case  of  equation  (2-16)  the  relation  is 


r   T  —MIX  T       -MTY    1 

iny        =  {(L  AG   )   -  (L  lira  AG^^^)  )/RT  -  S,nx     (2-22) 
oa a    -  o  —    a         oa 

-(a) 


In  terms  of  the  equilibrium  constants  of  a  single  set  of  reactions. 


W^G^  =  G^  +  RT  5^5,nK  +  RT(£nx  -  l/^-nx)  (2-23a) 


or  in  terms  of  activity  coefficients 


W  £nY  =  £nY  +  E,   inK  +  inx     -  w'^Jlnx  (2-23b) 

o  o 


Finally,  for  the  entropy. 


T-E   -T7       T  ^^^^ 
-o     -   9T 


+  R(ilnx  -  W^£nx)  (2-23c) 


P,N 
~o 


o 


T 
Of  course,  for  successive  reactions,  ^  JlnK  must  be  replaced  with  the 

general  expression. 

Unfortunately,  due  to  the  presence  of  the  logarithm,  the  excess 

function  expressions  do  not  have  the  kind  of  cancellation  there  is  for 

the  partial  molar  and  mixing  properties. 
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An  apparent  simplification  of  equations  (2-20)  is  normally  made 
for  electrolyte  solutions.   A  set  of  components  is  defined  for  which 


T 
In  X  =  W  In  X  (2-24) 
h    — n  


T      T 
where  W  is  W  with  its  elements  suitably  normalized, 
n 


Then, 


T 
In  Y^  =  W  In  Y  (2-25) 


However,  this  results  in  a  set  of  components  which  are  not  those 
which  could  be  used  to  make  up  the  system.   That  is. 


n 
o 


x_^   7*  X   E  N   /  y  N  ^ 
±a         oa    oa  qZ->    of 


This  difference  requires  extra  care  in  expressing  the  properties  of 
partially  dissociating  systems  (weak  electrolytes) . 

Complete  Dissociation 
An  important  special  case  of  the  above  formulae,  applying 
especially  to  strong  electrolytes  and  "solutions  of  groups"  methods, 
is  the  case  of  complete  dissociation.   Here,  we  simplify  W  by  elimi- 
nating rows  of  zeros,  making  it  into  a  matrix  V.   That  is,  the  reactions 
are  of  the  form 


B  >  Z(v).  .  P. 
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The  V  matrix  contains  the  stoichiometric  coefficients  of  only  the 
products,  and  subscripts  for  reactions  have  the  same  index  as  those 
for  the  corresponding  components.   Solutes  or  non-dissociating  compo- 
nents are  considered  to  undergo  the  reaction 


B   >  B. 

J  J 


so  that  they  appear  in  the  final  expressions. 
Our  equilibrium  equations  now  become 


y  =  v^y  (2-26a) 

o 


M  =  V^M  (2-26b) 


and  the  material  balance  is  written 


N  =  V  N  (2-27) 

—      o 


Now,  however,  the  species  vectors  do  not  include  any  dissociating 
components,  so  are  of  dimension  n  -  n  +  n  ,  where  n   is  the  number 

OSS 

of  components  that  do  not  undergo  dissociation  (solvents  =  n   <  n  ) . 
In  the  case  of  ions  and  groups,  the  actual  standard  states  are 
left  to  be  defined,  so  we  choose  JJ   as  for  equations  (2-15) 


y°  =  vV  (2-28) 

o 
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with  no  loss  of  generality.   This  leads  to  a  projectable  vector  of 
activities. 


v^^na  =  £na  (2-29a) 


or  V  Ir^     =  £nY  +  £nx  -  v^lnx  (2-29b) 

o      o 


As  noted  above,  the  "mean"  component  set  can  be  chosen   equation  (2-24) 
so  that  for  completely  dissociated  salts  equation  (2-25)  becomes 


ilny  =  V  ^^  (2-30) 


with  mean  ionic  activity  coefficients  Y,  and  mean  ionic  concentra- 
tions for  x+  (Harned  and  Owen,  1958;  Robinson  and  Stokes,  1965). 

Application  to  Analytic  Solution  of  Groups:  UNIFAC  and  UNIQUAC 

There  are  two  general  situations  where  group  contribution  methods 
are  used  in  correlation  of  thermodynamic  properties.   One  is  when  the 
parameters  of  a  molecular  correlation  are  obtained  by  summing  contribu- 
tions from  the  atomic  groupings  in  the  molecules.   For  example,  energy 
and  volume  parameters  have  been  calculated  in  various  models  from  sums 
of  group  energy  and  volume  contributions  (Reid  et  al. ,  1977).   The  method 
of  projection  operators  as  we  discuss  above  does  not  add  anything  to 
these  developments  because  the  groups  are  not  considered  to  have  any 
thermodynamic  property  values.   The  projection  of  parameters  does  not 
require  the  same  care  as  the  projection  of  properties. 
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The  other  situation  where  group  contributions  are  used  is  in  such 
cases  as  the  "analytical  solution  of  groups"  methods  (ASOG)  (Derr  and 
Deal,  1968,  1969,  1973).   Here  the  groups  are  considered  to  have  activity 
coefficients  which  are  summed  to  obtain  molecular  activity  coefficients. 
The  present  analysis  leads  to  significant  results  because  it  allows  a 
detailed  exploration  into  the  relationships  between  such  theories  for 
groups  (species)  as  UNIFAC  and  their  complementary  molecular  (component) 
theories  such  as  UNIQUAC.   Projection  operators  show  how  it  is  that  the 
two  equations  do  not  yield  the  same  composition  dependence  for  system 
partial  molar  properties.   As  shown  below,  this  difference  can  yield 
complete  miscibility  using  UNIFAC  but  only  partial  miscibility  for 
UNIQUAC  when  the  two  binary  parameters  for  the  latter  are  obtained  only 
from  the  infinite  dilution  activity  coefficients  predicted  by  UNIFAC. 

The  UNIQUAC  model  (Abrams  and  Prausnitz,  1975;  Maurer  and  Prausnitz, 
1978)  for  mixtures  gives  an  expression  for  the  excess  Gibbs  free  energy 
of  a  system  of  molecules.   In  the  UNIQUAC  Functional-group  Activity 
Coefficients  model  (Fredenslund  et  al. ,  1975,  1977a, b;  Skjold-Jorgenson, 
et  al. ,  1979)  (UNIFAC),  the  residual  portion  of  the  excess  Gibbs  free 
energy  is  written  in  terms  of  the  group  contributions  using  the  identi- 
cal expressions  for  the  composition  relations,  but  involving  projections 
of  the  group  parameters. 

This  expression  was  chosen  by  physical  argumients  after  UNIQUAC  was 
developed  and  the  question  can  be  raised  about  whether  equation  (C-3b) 
represents  a  proper  projection  of  (C-4) ,  which  would  then  lead  to  iden- 
tical results  for  both  equations.   If  not^  the  activity  coefficients 
can  be  matched  at  only  one  composition  other  than  pure  component. 
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Since  it  is  preferable  to  project  a  complete  thermodynamic  property 
instead  of  a  portion,  we  have  examined  the  energies,  instead  of  the  free 
energies,  as  given  in  the  derivations  of  the  models.   This  eliminates 
the  combinatorial  term  which  yields  the  same  results  for  both  equations 
and  focuses  only  on  the  residual  terms  where  the  differences  arise. 

Figure  4  shows  how  projections  can  be  applied  to  solutions  of 
groups.  A  three  dimensional  representation  is  given  of  the  excess 
enthalpies,  H  ,  of  a  solution  of  hydroxyl,  methyl  and  methylene  groups 


E    E       r 

H   :^U   EU-Tx. 
.   1 


lim  U 

x.^1  (2-31) 

1 


where  equation  (C-6)  has  been  used  with  parameters  AU    from  Skjold- 
Jorgensen  et  al.  (1979).   Vertical  planes  representing  mixtures  of 
n-pentane  with  methanol,  ethanol  or  n-pentanol  are  shown  with  the 
projections  being  the  heavy  lines  in  the  curved  surface.   The  ideal 
solution  line  for  these  binary  mixtures  is  the  straight  line  in  the 
plane  which  connects  the  values  of  H   at  the  two  pure  component  posi- 
tions.  A  portion  of  the  methanol-n-pentane  line  is  shown.   The  desired 

value  for  the  system,  H  ,  is  the  vertical  distance  from  the  ideal  solu- 

E  E 

tion  line  to  the  H   surface.   As  can  be  seen,  the  shapes  of  the  H 

curves  will  vary  because  curvature  of  the  H   surface  varies  with  chain 

length. 

While  the  plot  represents  an  exact  projection,  the  ASOG  method 

also  relates  the  group  and  molecular  properties.   In  fact,  because  the 

same  analytic  expressions  are  used,  expressions  for  the  partial  molar 

energies,    U        and  U„ ,    look  very   similar, 
ori  ii 
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EXCESS    ENTHALPIES    FROM   UNIFAC    FOR   n-ALCOHOL-n-PENTANE    SYSTEMS   AT    298K 

n-PENTANOL 


n-FENTANE 


CH. 


Figure  4.      Excess   enthalpies    from   U?1TFAC    for   n-nlcohol-n-pcntnne 
systems   at    298K. 
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However,  the  test  can  be  made  whether  the  true  projection 


U   =  I  V.  U„  (2-32) 

on  '^    l^  a 


exists  as  required  by  equation  (2-26b).   The  actual  expressions  are 
given  in  Appendix  C.   They  show  that  due-  to  the  complicated  interrela- 
tions of  the  indices  in  what  are  otherwise  very  similar  summations,  it 
is  clear  the  molecular  energy  parameters  cannot  be  related  to  the  group 
energy  parameters  in  order  to  make  the  two  expressions  become  equivalent, 
This  is  because  the  weighted  group  surface  fractions,  0^.,  cannot  be 
projected  to  give  a  set  of  weighted  molecular  surface  fractions,  9   , 
consistent  with  the  projections  of  the  group  surface  areas,  Q,  to  the 

molecular  surface  areas  q,  and  N   to  N.   It  is  possible  to  achieve 

— o    ~ 

numerical  equivalence  of  the  energies  at  any  particular  solution  compo- 
sition besides  unity,  but  the  parameter  relationships  will  vary  with 
concentration.   This  means  that  the  composition  dependence  of  the  com- 
ponent activity  coefficients  calculated  from  UNIFAC  will  be  different 
from  those  calculated  from  UNIQUAC  even  if  the  values  match  at  any  one 
point,  say  infinite  dilution.   Thus,  UNIFAC  should  be  considered  as 
different  from  UNIQUAC  as,  for  example,  the  Wilson  equation  is. 

If  the  numerical  differences  are  not  large,  then  the  procedure 
suggested  by  Fredenslund  et  al.  (1977b)  can  be  used  to  minimize  computa- 
tion time  by  finding  UNIQUAC  component  binary  parameters  from  UNIFAC 
infinite  dilution  activity  coefficients  and  calculating  vapor-liquid 
equilibrium  using  UNIQUAC  equations. 
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On  the  other  hand,  a  particularly  significant  problem  can  occur 
if  the  difference  is  large.   Since  the  slope  of  the  binary  UNIQUAC 
activity  coefficient  is  always  less  than  that  predicted  by  UNIFAC  at 
low  concentrations,  when  the  infinite  dilution  values  are  matched, 
situations  can  arise  where  UNIQUAC  predicts  liquid-liquid  immiscibility 
even  though  UNIFAC  does  not.    For  instance,  this  happens  in  the 
isopropanol-water  system  at  1.013  bar  (see  Figure  5).   This  problem 
can  be  avoided  by  fitting  the  UNIFAC  infinite  dilution  values  to  obtain 
Wilson  equation  parameters.   In  fact,  this  may  be  the  best  procedure 
in  any  case  for  the  miscible  systems. 


Multiple  solutions  for  UNIQUAC  parameters  can  be  found  by  this 
procedure  in  some  systems  (Fredenslund ,  1980),  but  we  have  not  found 
this  in  the  present  case. 
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Figure  5.   Free  energies  of  mixing  for  isopropanol(l)-water (2)  at  1  atm. 


CHAPTER  3 
STATISTICAL  THERMODYNAMICS  OF  "REACTIVE" 
COMPONENTS  IN  THE  GRAND  ENSEMBLE 


Introduction 
There  are  two  general  methods  for  relating  thermodynamic  properties 
of  solutions  to  statistical  mechanical  correlation  functions.   One  arises 
from  the  canonical  ensemble,  usually  assuming  pairwise  additivity  of  the 
intermolecular  potentials,  which  expresses  excess  energies  or  free  ener- 
gies in  terms  of  differences  of  integrals  involving  the  species  pair 
potentials  and  radial  distribution  functions  (Hill,  1956;  Reed  and 
Gubbins,  1973).   The  other  method,  from  the  grand  ensemble,  relates  con- 
centration derivatives  of  the  chemical  potential  to  integrals  over  the 
total  correlation  function  (Kirkwood  and  Buff,  1951;  Hall,  1971; 
O'Connell,  1971)  and  the  direct  correlation  function  (O'Connell,  1971). 
While  modelling  studies  have  been  done  for  both  procedures,  the  latter 
has  recently  been  shown  to  be  a  powerful  technique  for  solutions  con- 
taining supercritical  compounds  dissolved  in  liquids  at  elevated  pres- 
sures (Mathias  and  O'Connell,  1979,  1980a, b).   The  correlation  of  elec- 
trolyte solution  properties  (DeGance  and  O'Connell,  1975)  is  another 
application  of  interest  and  will  be  treated  in  subsequent  chapters. 
Furthermore,  the  direct  correlation  functions  from  the  grand  ensemble 
can  be  used  in  the  formalism  of  Chapter  2.   The  purpose  of  this  chapter 
is  to  apply  the  projection  formalism  in  such  a  way  as  to  establish 
fundamental  statistical  mechanical  relations  in  systems  with  arbitrary 

reaction  schemes. 
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Components  and  Species  in  Complex  Systems 
We  consider  the  closed  system  of  Chapter  1,  constructed  from  n 
neutral  components  that  may  undergo  reactions,  dissociations  and  asso- 
ciations to  arrive  at  an  equilibrium  system  of  n  chemically  identifiable 
or  hypothetical  species,   neutral  or  ionic.   As  before,  the  material 
balance  which  relates  the  vector  of  the  particle  numbers  of  the  species, 
N,   to  that  of  the  components,  N  ,is 


N  =  W  N  (3-1) 

—   —  — o 


where  W  is  as  defined  previously. 

In  this  more  detailed  analysis,  we  must  use  the  structure  of  W 
and  related  matrices  to  obtain  the  desired  results.   We  look  here  at 
the  case  of  no  successive  reactions,  and  leave  the  general  case  to 
Appendix  D.   In  this  simple  case,  we  have,  as  in  Chapter  1, 

W  =  L  +  V  5  (3-2) 


1  In  this  chapter,  the  vector  N  will  have  two  meanings.   When  used 
in  an  obviously  macroscopic  sense,  N  is  the  vector  of  the  particle  num- 
bers of  the  species  in  the  system  of  interest.   However,  in  the  statis- 
tical mechanical  derivations  beginning  in  the  next  section,  N^  will  de- 
note the  composition  vector  of  a  specific  system  in  the  ensemble,  and 
<N>will  denote  the  average  of  N  over  all  the  systems  of  the  ensemble, 
and  so  corresponds  to  the  N  used  previously.   Later,  we  will  return  to 
using  N  instead  of  <N>,  which  is  cumbersome. 
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The  L  matrix  is  of  dimension  n  x  n  ,  where  the  first  n   rows  have  the 
~  o  o 

form  of  an  n  x  n  identity  matrix  and  the  remaining  n-n  rows  contain 
o    o  o 

only  zeroes.   L  can  then  be  written  as 


L  = 


0 


(3-3) 


The  _v  matrix  is  of  dimension  n  x  r  and  can  be  decomposed  in  the 

following  manner.   We  call  the  first  n  rows  v^  (n  x  r)  which  contain 

o      —Co 

the  stoichiometric  coefficients  of  the  components  in  the  reaction 
scheme.   We  next  note  that  the  reaction  equilibria  are  defined  by 

v\  =  0  (3-4) 

where  v_   again  is  the  vector  of  species  chemical  potentials.   Since  the 

n  component  chemical  potentials  are  independent  and  equation  (3-4)  puts 

r  contraints  on  p ,  a  given  system  (fixed  p  )  has  n  +r  total  contraints 
—  — o       o 

on  p.   This  leaves  n-n  -r  "degrees  of  freedom"  in  y.   That  is,  n-n  -r 
—  o  —  o 

chemical  potentials  can  be  varied  independently  without  changing  the 
values  of  P   and  without  violating  equation  (3-4).   Of  course,  the 

other  r  chemical  potentials  will  be  functions  of  the  chemical  potentials 

2 
which  are  chosen  to  be  "independent". 

Following  these  identifications,  we  now  define  matrix  \)   (r  x  r) 
which  contains  the  stoichiometric  coefficients  of  the  "dependent" 
species,  and  the  matrix  v  which  contains  the  stoichiometric 


^These  are  not  degrees  of  freedom  in  a  normal  sense  in  that  a 
fixed  T,V,Pq  fixes  all  of  \i_,   but  if  jj_q  is  fixed  there  may  be  multiple 
solutions  to  v^\i_   =  0. 

^The  identification  of  these  species  is  not  unique.   There  will  be 
some  species  which  cannot  be  Independent  in  this  sense,  but  in  general 
there  are  many  ways  to  divide  the  independent  and  dependent  species. 
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coefficients  of  the  "independent"  species.   Because  of  the  manner  of 
construction,  v  will  be  nonsingular,  so  v   will  always  exist. 
Using  these  identifications,  we  can  rewrite  v  and  W  as 


^D 


^I 


(3-5) 


and 


W  = 


(3-6) 


The  matrix  W  has  dimensions  n  x  n  with  rank  n  ,  provided  all  the 
-  o  o 

original  components  appear  in  the  final  system.    This  means  that  W  maps 

N  onto  an  n  -dimensional  "component  subspace"  in  the  n-dimensional 
-o         o 

"species  space,"  and  therefore  there  exists  a  space  of  dimension  n-n 

in  the  species  space  which  is  orthogonal  to  the  component  subspace,  and 

thus  to  the  columns  of  W.   We  define  Z  as  an  n  x(n-n  )  matrix  made  up 

-  -  o 

of  a  set  of  vectors  which  span  this  space  orthogonal  to  W.   Then 


Z^  =  0  and  so,  Z  N  =  0 


(3-7) 


Of  course,  any  set  of  vectors  which  span  this  space  can  be  used  to 
construct  Z,  and  in  general  each  element  of  Z  is  a  function  of  the 


1 Another  situation  arises  in  the  case  of  completely  dissociated 
electrolytes,  for  example.   Then  the  rank  of  W  can  be  less  than  n^. 
The  formulation  must  be  slightly  modified  so  that  the  above  analysis  is 
still  applicable.   The  details  are  given  below. 


37 


extent  of  reaction.   For  convenience,  we  will  construct  a  specific  Z 
consisting  of  two  parts:   Z  ,  which  is  n  x(n-n  -r)  and  has  elements 
which  are  independent  of  the  C  's,  and  Z  which  is  n  x  r  and  has 

elements  which  are  functions  of  the  5,  's.   Z  can  be  written  explicitly 

k     -o 

as 


Z  = 
-o 


0 


,  T,-l   T 


(3-8) 


wit 


th  0  being  n  x(n-n  -r)  and  I  being  (n-n  -r)  x  (n-n  -r) .   It  is  clear 


o      o 


that 


Z^W  =  0^  (I  +  v^p  -  VjV'^D  +  IVj  =  -Vj  +  V;,  =  0 


as  desired. 


The  Z  matrix  is  defined  by  the  following  relations: 


Z  W  =  0 


(3-9a) 


Z,  Z  =  0 
-1  -o 


(3-9b) 


5l  ?1  =  i 


(3-9c) 


which  completely  define  Z  .   The  Z  matrix  can  now  be  written  as 


Z  =  (Z  ,  Z  ) 
-    -o  -1 


(3-10) 


which  does  span  the  (n-n  )-dimensional  subspace  orthogonal  to  W  as 

'^  o 

desired. 
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We  also  must  construct  two  n  x  n  matrices,  U  and  Y,  which  project 
in  the  direction  opposite  to  W.   First,  U  is  important  since  it  spans 
the  same  space  as  W,  whereas  Y  is  important  because  its  elements  are 
independent  of  the  extents  of  reaction.   The  defining  relations  for  U  are 


T 
U  W  =  I 


(3-lla) 


T 
U  Z  =  0 


(3-llb) 


and  it  has  the  property  that 


u\  =  N 
-  -   -o 


(3-12) 


We  can  write  Y  explicitly  as 


0 


(3-13) 


with  I  being  n  x  n  and  0  being  (n-n  -r)  x  n  .   Also 
—        o    o     -  o       o 


T 

Y  N  =  N 
-  -   -o 


(3-14a) 


but 


Y^Z  +   0 


(3-14b) 
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Grand  Ensemble  for  Systems  of  "Reactive"  Components 
To  establish  the  statistical  thermodynamic  relations  among  the 
components  and  species,  we  begin  by  establishing  an  ensemble  with  /7_ 
systems,  total  energy  N  <E> ,   and  /V  <N.>  particles  of  each  species  i. 

The  number  of  ways,  t,  of  having  a  given  distribution  of  m   's,  where 

Nq 

n        is  the  number  of  systems  in  quantum  state  q  with  particle  number 
vector  N  is 


n  n  n  ! 

q  N   -^ 


In  Stirling's  approximation, 


£nt  =  N^ZnN^  -  N^  -   H    (n   £nn^  -  n   )  (3-16) 

N  q 


Using  the  assumption  of  equal  a  priori  probabilities,  int  is  maximized 
to  find  the  most  probable  distribution,  subject  to  the  constraints  of 
constant  total  size  and  energy  and  the  overall  material  balance 
equations.   These  constraints,  with  their  Lagrange  multipliers,  are 


lln       =N^  a  (3-17a) 

N  q   - 


N  q   -   -^ 
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y  y  n   U  N  =  /V  U  <N>  =  Nn              ln\  (3-18a) 

^  '^  Nq  -  -    T-   -     T-o     


y  y  n  Z^^  =   N  z'^<N>  =  0        2.ne       (3-18b) 

Z,  Z,   Nq T-   -    -        

N  q   - 


where  toX  and  toe  are  vectors  of  dimension  n  and  n-n^  respectively  and 

whose  elements  are  ZnX    and  Jinev  The  equation  for  the  maximum  in  tot  is 

a  K. 

d(iLnt)   +  a   I  Idn       -   &   I   Ie^     dn       +   (^A)^   I   1"^% 
N   q      -  N   q   -^        -  N   q 


+   ilnO)'^   I   Izhdn^     =   0  (3-19) 

N    q 


Equation    (3-16)    gives,   upon  differentiation 


d(£nt)    =  -   y   y    Inn       dn 
^  t,  Nq        Nq 

N   q  -  - 


which  with  equation  (17)  leads  to 


-inn,,     +  a  -   3E,,  +  (£nA)VN  +  (!lne)^z'^N  =  0       (3-20) 
Nq         Nq    -  -    — ^   -  - 


since  in  the  method  of  Lagrange  multipliers  the  n       are  treated  as 
independent. 

Solving  equation  (3-20)  for  n      ,   we  have,  for  all  values  of  N  and  q 


n        =   exp  {a  -  3E„  +  (£nX)VN  +  (to0)'^z\}        (3-21) 
Nq  Nq    -  -    -  - 
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Defining 

5(A_,e^,6,V)  =  I   I  ex-p   {-6E^  +  (£nA)VN  +  (£ne)^z'^N}     (3-22) 

the  probability  of  finding  a  system  with  particle  number  vector  N  in 
state  q  is  written  as 


^q   exp  {-BEjj  +  (£nA)V  +  (£n9)^Z^  N} 

PNq  =  ^-  =  = = =— ^       (3-23) 

-^    T  E(X,9,B,v) 


In  the  usual  case  of  nonreacting  components,  the  partition  function 
would  be  written 


:(y,3,V)  =   II   exp   {-3E   +  3y  N}  (3-24) 

N  q         -^ 


where  N  is  the  same  vector  as  in  equation  (3-22)  ,  3  is  1/k  T  and  p  is 
—  B      — 

the  vector  of  species  chemical  potentials.   By  comparing  equations  (3-22) 
and  (3-24)  we  should  expect  that  £nA  and  fcn6  can  be  identified  by  the 
relation 

3y  =  U  _2riA  +  Z  _Kjie  (3-25) 

This  identification,  if  correct,  has  several  interesting  results. 
Using  equation  (2-5)  and  equations  (3-11) 


3y  =  3w'^y  =  JlnA  (3-26) 

-o    -  -   
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Furthermore,  equations  (2-2,  (2-3)  and  the  Gibbs-Duhem  equation 
lead  to,  after  some  manipulation  and  use  of  equations  (3-11) 


{(£nX)'^dU^  +  (S,ne)^dz'^}  <N>  =  0 


(3-27) 


provided  the  differentials  are  taken  along  an  equilibrium  path. 

Equation  (3-26)  identifies  B   £nA  with  the  chemical  potentials  pf 
the  components.   Equation  (3-27)  will  be  useful  in  later  manipulations. 

The  identification  in  equation  (3-25)  is  only  a  proposition,  and 
it  must  be  shown  to  be  a  solution.   We  will  assume  it  is  a  solution  and 
then  test  the  validity  of  this  assumption.   We  begin  by  examining  the 
average  entropy  for  a  system  in  the  grand  ensemble,  ■ 


cS>  =  -k   y  y  p   2,np., 
^  N  q   -"^   -"^ 


or,  substituting  equation  (3-23)  for  p^,  in  this  relation 

Wq 


<S>/,   =  6<E>  -  (£nA)  N  +  Znl 
kg  -o 


(3-28) 


For  the  Helmholtz  free  energy.  A,  we  have 


BA  =  B<E>  -  <S>/,   =  (£nA)  N  +  inl 
kg    -o 


(3-29) 


The  chemical  potential  of  component  a  in  component  space  is 


By_.  =  B- 


oa 


aA 

3N 


oa 


9(£nA)' 


3N 


T,V,N 


oa 


N     +   inX 
-o  a 


9£nl 


9N 


oY?^a 


T,V,N 


oa 


(3-30) 


oy^a 


T,V,N 


oy/ot 
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where  it  is  important  to  note  that  it  is  the  amounts  of  the  other  com- 
ponents (N  )  as  initially  charged  to  the  system  before  reaction  occurs, 
and  not  the  amounts  of  the  species,  that  are  held  constant.   We  assume 
that  a  change  in  a  component  concentration  is  done  with  all  reactions  at 
equilibrium.   Thus  the  species  concentrations  are  known  from  the  com- 
ponent values,  N  ,  and  the  material  balance  based  on  equilibrium  extents 
of  reaction   (W) .   With  these  restrictions  all  the  derivatives  in  equa- 
tion (3-30)  are  taken  along  an  equilibrium  path,  making  equation  (3-27) 
applicable. 

Since  £.nA  and  £n9  are  functions  of  N  ,  equation  (3-30)  can  be 
written 


d(in\)' 


PP, 


Oct 


3N 


oa 


8(5,nA) 


N  +  InX  ^- 
-o        01 

T,V,N 


9N 


oa 


y  <N>  -  (lnA)^— 
T,V,N 


Y?^oi 


oa 


<N> 
T,V,N, 


3   , 


3(£ne) 


3N 


oa 


T  T  ^5 

Z  <N>  -  (£nA) 


3N 


T,V,N 


oa 


<N> 
T,V,N^ 


3   / 

yra 


(3-31) 


The  first  and  third  terms  on  the  right  hand  side  of  equation  (3-31) 
cancel  by  equation  (3-12) ,  the  fourth  and  sixth  are  zero  by  equation 
(3-27)  and  the  fifth  term  is  zero  by  equation  (3-7).   This  leaves  us 
with  the  desired  result 


oa      a 


or 


3y  =  ^nA 
-o   


(3-32) 
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as  proposed  in  equation  (3-25)  and  (3-26).   Thus  we  find  that  these 
equations  are  indeed  a  solution  for  the  Lagrange  multipliers  HnX   and 

£n9. 

Equation  (3-32)  not  only  defines  the  Lagrange  multiplier  2,nX 

a. 

but  also  shows  us  what  the  independent  variables  for  the  grand  canon- 
ical ensemble  are.   We  already  know  that  a  macroscopic  (canonical)  sys- 
tem can  be  deteinnined  by  fixing  T.V  and  N  ,  a  total  of  n  +2  variables. 
Similarly,  in  the  grand  ensemble,  we  fix  the  same  number  of  variables 
T,  V  and  y   (or  £nA)  to  completely  determine  the  thermodynamic  state. 
This  fact  is  very  important  in  taking  derivatives  of  ensemble  averages, 
since  this  means  there  is  a  relationship  between  the  elements  of  g,n6 
and  the  elements  of  ZnX.      This  dependence  and  important  thermodynamic 
derivatives  of  the  ensemble  are  given  in  Appendix  E. 
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Composition  Fluctuations 
Macroscopic  thermodynamic  properties  of  the  ensemble  can  be  re- 
lated to  microscopic  properties  through  the  density  fluctuations  in  the 
systems  of  the  ensemble.   We  define  the  fluctuation  matrix.  A,  by 


with 


A  =  -^   (<NN^>  -  <N><n'^>) 


(3-33) 


1  T 

<N>  =  -=-  I  I  N  exp  {-3Ej^  +  N  (U  £nX  +  Z  £n6)}  (3-34a) 

"   N  q  - 


<N>  =  y  <N.> 


(3-34b) 


<NN  >=-r  5!  I  NN   exp  (-3E   +  N  (U  laX  +   Z  £n9)  }   (3-34c) 
"  N  q  - 


Equation  (3-34a)  can  be  differentiated  with  respect  to  a  component 
chemical  potential  along  the  equilibrium  path  to  obtain  the  relationship 
between  the  fluctuations  and  macroscopic  properties  (details  are  given 
in  Appendix  E) .   This  differentiation  leads  to 


1  8<N> 


<N>93m 


oa 


1    T 
=  -=■  <NN  > 
<N>  — 


T,V,ii 


oY?*a 


8Mj 

(Y)   +  Z--i 

-  a   -ody 


oa 


T,V,y 


OY?^ot 


■^  <N><N^>  (Y) 
<N>  -  -    -  a 


(3-35a) 
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=  A   (Y)   +  Z  ^ 
-  ct   -o  9y 


(3-35b) 


oa 


T,V,u 


oY?^a 


where  the  notation  (Y)   denotes  the  a-th  column  of  Y. 

-  a 

We  can  use  equation  (3-14a)  to  project  the  left  hand  side  of  this 
equation  into  component  space.   Defining  the  matrix  A 


9N 


<-i\i 


-■^ga  "  <N>  93m 


oB 


oa 


(Y) 


T  9<N> 


<N>  -'3  93y 


T,V,u 


oa 


oy/a 


T,V,y 


oy^a 


1     T    T 

<i  (^3  ^^  =• 


9jjj 


(Y)   +  Z   , 

-  a   -o  9y 


oa 


T,V,y 


OY?^a 


<i  (Dg  <n><n''>  (y)^ 


(3-36) 


Using  equation  (3-7)  we  can  rewrite  ^  as  a  projection  of  A, 


^-\,  -   (^->3  '- 


(Y)   +  Z 


9jjj 


a   -o  9y 


oa 


T,V,y 


oy^a^ 


T    ^^1 


-o 

T,V,y 


oy^a 


A  (Y) 

-  a 


(3-37) 


The  next  section  shows  how  the  A  matrix  is  related  to  the 
microscopic  properties  of  the  species,  thus  forming  a  bridge  between 
them  and  the  macroscopic  component  properties. 
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Correlation  Function  Integrals 
The  method  of  Kirkwood  and  Buff  (1951)  can  be  used  to  relate  the 
fluctuation  properties  of  integrals  of  the  radial  distribution  function 
and  direct  correlation  function.   We  define  spacially  and  orientation- 
ally  dependent  one  and  two  particle  densities  by 


N, 


:(1) 


f,  {v-,;(xi-.)    =  y   (r  -r-,)  6(a)  -w, )    i,  =  l,...,n    (3-38a) 
1   -1  -1     ^-1   -p  -1    -p  -1 


p=l 


:(2) 


N.   N. 

1     X 


fl-^£l,£2'^l'^2)  =    \      \    '%-'-l^    '%-'-2^    '%-^-l^ 
-J  p=l  q=l 


p=l  q= 
(p^q  if  i=j) 


X  6(cj  -a)„) 
-q  -2 


(3-38b) 


where  6(r-r^)  is  the  three  dimensional  Dirac  delta  function,  r  is  the 
position  vector  and  o)  is  the  vector  of  Euler  angles  of  molecular  orien- 
tation. 

By  definition,  we  have 


Vff^(£r-V  dr,d.^  =  N. 


(3-39a) 


:(2) 


V^ii   (£l'l2'^l'^2>  'i£l'^£2^-ld-2  =  Vj  -  ^i^ij  ('-^'"^ 


We  now  define  ensemble  average  number  densities  by 


pf  ^(Hr^i>  '  <4'^(£r-i» 


(3-40a) 
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(2) 


(2), 


Pl^^^l'll'^^V^l^   =  <fli  ^ivii'^v'^i^^ 


(3-40b) 


Equations  (3-39)  are  integrated  to  obtain 


V 


(1) 


p)^\r-,i>i,)    dr.du).  =  <N,> 


(3-41a) 


Q 


(2) 


Pi-  ^^l'£2'-l'-2^  '^-l'^-2^-l^-2  "  ^^i^j^  "  ^^i'^'^ij    (3-41b) 


For  fluid  systems,  we  define  the  mean  density  and  the  pair  correlation 
function  by 


(1), 


p^  =    fip^  (li'^'^i)    =   <N^>/V 


(3-42) 


and 


2  (2) 


PiPi^i1^-12'-l'-2^  =  ^  P^j  (£i'£2'-l'-2^ 


(3-43) 


where  f2 


dco. 


The  pair  correlation  function  can  be  integrated  over  its  arguments 

and  related  to  the  concentration  fluctuations  using  equation  (3-3%),  as 


p  p 

<N.N.>  -  <N.>6.  .  =  -^^   V 


§•  •  (£l2'-l'-2^  d£j^2'^-l'^-2    (3-44a) 
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or 


<N.N.>  -  <N.>6.  .  =  p.p .V 
1  J      1   13     12) 


/gij^^l2'^l'^2^\^f7/Hi2 


(3-44) 


The  angle-averaged  (or  "centers")  total  correlation  function,  h.,(r^„), 
is  defined  by 


(3-45) 


Then  using  matrix  notation 


(H)y 


<N> 
V 


Nj(^12^'^Hi2 


=  <N> 


<N.N.> 
1  J 


<N.><N.> 
1    J 


<N.><N.> 
1    J 


6.. 
13 

<N.> 

1 

(3-46) 


we  f  ind 


<N.><N.> 

<N  >       -^ 


<N.N.>  -  <N.><N.>    <N.> 
J-  3 1   3 i_ 


<N> 


<N>   ij 


(3-47) 


letting  X.  =  <N.>/<N>,  equations  (3-33)  and  (3-47)  yield 


(A)..  =  X.6..  +  x.x,(H).. 
^-'13    1  13    1  3  -  13 


or,  defining  a  matrix  (X) . .  =  x.S. . 
^  -  13    1  13 


A  =  X(I  +  H  X) 


(3-48) 
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A  "centers"  direct  correlation  function  can  be  defined  by  the 


equation 


n 


^ilj^^l3^PAii^^23)'^H3 


(3-49) 


Integrating  over  r   and  multiplying  by  <N>/V  we  find 


<N> 


h.  .dr 


<N> 


V 


ij  -12    V   j 


<N> 


n 


c.  .dr  T  +  ,, 
V  ^J  -^2    ^   £=1  ^V 


^''ii^^-U^f'i^l^^-li^ 


X  dr3dr^2 


(3-50) 


For  systems  where  the  matrices  A,  H  and  C  are  all  bounded  and  non- 
singular  ,  Fourier  transform  theory  (Pearson  and  Rushbrooke,  1957)  can 
be  used  to  obtain 


H  =  C  +  HXC  =  (I+HX)  C 


(3-51) 


where 


(C)..  = 


<N> 


-'ij    V 


c.  .  (r  ^)dr 
,  ij  -12  -12 


Some  manipulation  of  equations  (3-48)  and  (3-51)  leads  to  the  relation- 
ship between  C  and  A, 


iHere  it  has  tacitly  been  assumed  that  the  integral  of  Cij(r;p)  is 
not  divergent.   This  is  not  the  case  for  certain  systems,  such  as  those 
with  coulombic  interactions.   In  such  a  case,  the  derivation  of 
Chapter  4  must  be  used. 
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A(X  ''"-0=1 


(3-52) 


Since  we  have  only  the  components  as  degrees  of  freedom,  we  want  to 
project  this  into  component  space.   Following  equation  (3-37)  we  pre- 
multiply  by 


T    ^Hi 
(Y)^  +  ^ 

—  a   9y 


oa 


-o 

T,V,u 


oy^^a 


and  post-multiply  by  W,  giving 


T    ^Hi 
(Y)'  +^ 

-  a   8y 


oa 


z" 

-o 


T,V,y 


oy/a 


A(X  ■*"  -  C)W 


^vl^- 


+ 


9Hi 


9y. 


oa 


T 
Z  W 
-o  - 


T,V,y 


oyt'cx 


(3-53a) 


or 


T 

3<N  > 


9  By 


oa 


T,V,y 


(X  -^  -  C)W  =  e 
-  -   -a 

oY^a 


(3-53b) 


where  e  is  a  vector  such  that  (e  )„  =  5   ,   Differentiation  of 

-a  -ct  p    aps 

equation  (3-1)  leads  to 


T 


96y 


oa 


T 
9N 
-o 


96y 


T,V,y 


oa 


T         ,T    9  W 
W     +  N 


-o    9  By 


oY^a 


T,V,y 


oa 


(3-54) 


oY?^ot 


T,V,y 


oY?^« 
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We  define   the  matrix  K     such   that 


o  9(W). 

(K   )  .      =      I     N        ■-      ^^ 

Y=l  oa 


(3-55) 


T,V,y 


OY?^ot 


Then  we  can  write  the  general  expression  in  component  space  by  com- 
bining  equations  (3-37),  (3-53b) ,  (3-54)  and  (3-55) 


^^■'■(X  ■''-C)W  +  K^(X  "'■-C)W  =  J 


(3-56) 


where  I   is  the  n  xn  identity  matrix.   Multiplying  through  by  4   ,  we 
have  the  desired  result 


4  "*■  =  (W  +  K)^(X  -"--OW 


(3-57) 


where 


o     9(W). 
(K).  =   5;   N   iX 


la    ^-,      OY   3N 
Y=l 


oa 


T,V,N 


oy/a 


n      n    . 

O        O     dp 


=   I   N  .  I 


oB 


^   oYo^T    9N 
Y=l    3=1     oct 


9(W).  I 


oY^3 


=  I      (A    '^)       (K   ).^   =    (A    V)  . 


(3-58) 


Equation  (3-57),  therefore,  relates  the  derivatives  of  the  component 
chemical  potentials  in  component  space  to  the  species  direct  correlation 


53 


functions  and  the  component  composition  derivatives  of  the  extents  of 
reactions  (through  W) . 

Because  of  its  complexity,  it  is  not  apparent  that  a  single  optimal 
procedure  exists  for  obtaining  changes  in  chemical  potentials  (and  also 
pressure)  by  integrating  equation  (3-57).   However,  in  the  case  of  to- 
tally dissociating  systems  this  is  straightforward. 
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Complete  Reactions 

One  kind  of  system  of  particular  Interest  Is  that  which,  for  each 
reaction,  there  is  a  limiting  reactant  which  appears  only  In  that  re- 
action and  whose  concentration  goes  to  zero  In  the  solution.    In  such 
a  case,  the  species  composition  vector  and  the  extents  of  reaction  can  be 
uniquely  determined  from  v  and  N  using  only  algebraic  relations  (no 
equilibrium  calculations).   Examples  of  this  kind  of  system  Include 
solutions  of  strong  electrolytes  and  "solutions  of  groups." 

Here  we  make  a  different  partitioning  of  the  matrices  of  interest. 
The  components  are  divided  into  two  groups:   the  limiting  reactants  which 
are  subscripted  L,  and  the  other  components  which  are  subscripted  S 
(solvents) .   The  species  which  are  not  found  as  initial  components  form 
only  one  group,  which  is  subscripted  P  (products).   The  W  matrix  is  then 


W  = 


I 

0   ^ 

"-S    ^ 

0 

I 

+ 

^-L 

^^-S    k> 

0 

0 

^-P 

(3-59) 


Since  a  system  of  only  complete  reactions  with  the  above  restrictions 
must  have  the  same  number  of  reactions  as  it  does  limiting  reactants,  we 
will  number  each  reaction  by  Its  corresponding  limiting  reactant.   If 
n^  is  the  number  of  solvents,  the  partitions  of  the  W  matrix  have  the 
following  dimensions: 


In  the  case  of  where  two  components  of  a  reaction  would  go  to 
zero  concentration,  one  is  chosen  as  the  limiting  reactant  and  the  other 
treated  as  an  ordinary  component.   This  situation  could  only  occur  if 
these  two  components  had  Initial  compositions  with  the  proper  stoichio- 
metric ratio.   In  the  derivation  given,  only  non-sequential  reaction 
schemes  are  considered.   This  is  general  because  any  sequence  of  com- 
plete reactions  can  be  written  as  a  non-sequential  scheme. 
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V  is  n  X  (n  -n  ) 

-S       S       OS 

v^  is  (n  -n  )  x  (n  -n  ) 
-L      OS      OS 

v„  is  (n-n  )  x  (n  -n  ) 
-P         O       OS 

^S  ^^  ^V^s^  ^  "s 


To  see  the  simplifications  that  result  in  the  case  of  complete 
reactions,  equation  (3-57)  can  be  written 


A^   =    (W+k'^)  (X  -"--O  (W+K) 


(3-60) 


This  equation  may  be  obtained  from  equation  (3-52)  by  following  a  pro- 
cess identical  to  equations  (3-53)  through  (3-58)  plus  postmultiplying 
equation  (3-52)  by  K  to  obtain 


T   -1        T    ^Hl 
^(W+K  ) (X   -C)K  =  Y  K  +  ^ 


-oa 


T 
Z  K  =  0 
-o- 


T,V 


(3-61) 


The  second  equality  comes  from 


and 


^Hi 


^Ho 


y'^K  =  0  and   Z^K  =  0 


is  the  Jacobian  matrix  of  the  derivatives 


9(Hx)i 


3y 


T.V 


oa 


T,V,y 


oY?^a 


Thus, 


(W+K)'^(X  """-OK  =  0 


(3-62) 
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Addition  of  equations  (3-58)  and  (3-62)  gives  equation  (3-60)  as  desired, 

The  projection  matrix  W  +  K  is  a  very  sparse  matrix  compared  to  either 

W  or  K. 

For  a  complete  reaction,  written 

n 

y   B.v.  =  0 
1=1 

with  component  a  being  the  limiting  reactant,  we  know  that 


N=0=N   +v£N  (3-63a) 

a       oa    aa  a  or 

a 

where  r  is  the  reference  component  for  reaction  a.   In  the  limit  of 
complete  reaction,  the  extent  of  reaction  can  be  written 

5  =  -  J-  ^  (3_63b) 

a     V   N 

aa  or 
a 


The  general  term  in  the  W  matrix  is 


n 

(^i)i3  =  ^i6+  ^°  >(-^N^)^Br  (3-64) 

a=n  +1     aa  or     a 
s  a 


Similarly,  for  the  K  matrix 


n      n  6  .N    -  N   6„ 

o      o     V.        ag  or     oa  Br 

(K).   =  I     N    I  (-- i^)6    °L « 

-  iB    ^   OY  ^    _,   V   '  yr        „2 
Y=l    a=n  +1   aa    a       N 

s  or 

a 


o   V.  o   V.   N 

ot^n  +1  aa       a=n  +i  aa  or 
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Addition  of  equations  (3-64)  and  (3-65)  yields  the  general  element  of 
the  projection  matrix 


o     V. 
la 


(W+K)   =  6   -  I  -J:2:  5 


a=n  +1  aa 
s 


(3-66) 


If  index  g  is  between  1  and  n   (solvent) ,  we  have  simply 

s 


(W+K) .„=6.„  l<B<n;l<a<n 

-   -  iB  i3  s' 


(3-67) 


which  defines  the  first  n  columns  of  (W+K) .   If  B  is  between  n  +1 

s  -  -  s 

and  n   (limiting  reactant)  we  have  three  possibilities:   For  i  being  a 

solvent  (1  5  i  <  n  ) , 

s 


(H+PiB=   V 


iB 


lliln;n+llB_n. 


BB 


s'  s 


(3-68a) 


For   i  being  a  limiting   reactant    (n  +1   <   i   <  n  ) 
°  so 


(W+K)^g  =  0 


n+l<i,B^n  i/ 

s  o 


(W+K)  .  „   =  1 —  =0  n+l<i,B<n  i  = 

-  -^2  ^B3 


(3-68b) 
(3-68c) 


Finally, for  i  being  a  product  (n  +1  <  i  <  n) 


(W  +  K) 


^f 


-'iB  V 


BB 


n+l<i<n;n+l<S<n  (3-69) 

o  s  o 


Since  v.     (made  up   solely  of   the  v      's)    is   diagonal,    equations    (3-65)    and 

L  pp 

(3-69)  can  be  written  in  matrix  form. 


W+K  = 


I 
0 
0 


-1 

—V  V 

-S-L 


-^P^L 


(3-70) 
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Partitioning  (X  -C)  in  a  similar  manner  yields 


x-i-c  = 


^x-i-c 
-s  -ss 


-^LS 


-£ps 


-LS 


-^LL 


-C 


PL 


-PS 

' 

-PL 

+ 

5'-2pp, 

V 

0 

0 

0 

h 

0 

0 

(3-71) 


■) 


The  first  matrix  on  the  right  hand  side  is  well-behaved  in  the  limit  of 
complete  reactions.   However,  the  second  matrix  diverges  and  must  be 
operated  on  by  (W+K)  before  the  limit  of  complete  reactions  is  taken. 
This  procedure  is  carried  out  in  Appendix  F  and  leads  to  the  equation 


f 


limit 
, complete  .   (W+K) 
reactions 


0 

0 

0 

0 

t 

0 

(W+K) 

=   0 

0 

0 

0 

(3-72) 


Equations  (3-60),  (3-70),  (3-71)  and  (3-72)  lead  to  the  following 
equation  for  A        in  the  complete  reaction  limit. 


-1 


5s'-^ss 


(^l')'^-4(5s'-?SS>  +  ^P^PP> 


^-(^s-Ss)^s-^^?s^p>^l' 
(^l'>'^4(^sSs)^s-^pSs^s 


^VgCj^Vp  +  vJ(x;^-Cpp)Vp}v-^ 


(3-73) 


Although  equation  (3-73)  appears  somewhat  complicated,  it  contains 
some  important  results  and  has  two  important  limiting  cases.   Most  im- 


portant of  these  is  that  there  are  no  contributions  to  A        from  direct 
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correlation  functions  involving  limiting  reactants.   Also,  it  can  be 
seen  that  A        is  symmetric  as  required. 

If  we  look  at  a  system  with  only  solvation  or  association  of  com- 
ponents, then  all  reactions  are  of  the  form 


n 


I   ^.T.B.   +  ^TtA.  -  B 


a=l 


aL  a     LL  L    P 


where  the  B  's  are  solvents,  B^  is  the  limiting  reactant  and  B_,  is  the 

a  LP 

solvated  product.   The v  matrix  is  the  identity  matrix  and 


A-^= 


^sHs 


^-(S'-^ss>^s  +  ^ls^\ 


-^s  Cps  +  (Xp  -Cpp)}vL 


(3-74) 


Equation  (3-74)  simplifies  still  further  for  systems  with  only  either 
solvation  or  association.   For  the  former,  v.  is  the  identity  and 


-1 


^s'-^ss 


-^^S'-^ss)  ^  5ps 


-(S'-Cs3)v3  +  4 


Vg^S^-^ss^^s  -  ^PS^S 


-VgCps  +  ^h'   -   Cpp) 


(3-75) 


For  the  latter  case,  v  =  0  and 


-1 


-s  2ss 


-1  T 


pT   -1 
^PS^L 


(v-')'(x;^pp)v-^ 


(3-76) 
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The  other  case  of  interest,  and  the  case  to  be  pursued  more  fully 

in  later  chapters,  is  that  of  complete  dissociations,  that  is,  reactions 

of  the  form 

n 

k=n  +1 
o 

where  B  is  the  limiting  reactant  (dissociating  component)  and  the 

Li 

B,  's  are  product  species.   The  solvents  do  not  take  part  in  the  re- 
actions (v„=0) .   In  this  case,  v  is  the  identity  matrix  and 


,-1 


-s  ^ss 

T 


T 
^PS^P 


(3-77) 


These  results  are  quite  simple  to  use,  in  fact  evaluation  of 
activities  requires  only  simply  quadratures,  not  the  solution  of  partial 
differential  equations  as  in  the  case  of  partial  reactions. 


CHAPTER  4 
SOLUTIONS  OF  STRONG  ELECTROLYTES 


Introduction 


O'Connell  and  DeGance  (1975)  give  a  formalism  for  expressing  the 
fluctuation  properties  of  electrolyte  solutions  (concentration  deriv- 
atives of  the  chemical  potential,  partial  molar  volumes  and  compressi- 
bility) in  terms  of  integrals  of  the  radial  distribution  functions 
(Kirkwood-Buf f  solution  theory,  Kirkwood  and  Buff,  1951;  O'Connell, 1971) , 
and  of  nondivergent  integrals  of  the  direct  correlation  function.   The 
formalism  was  based  on  a  constraint  on  these  radial  distribution  func- 
tions (Friedman  and  Ramanathan,  1970)  which  is  not  necessary  to  arrive 
at  the  desired  results  in  the  case  of  the  direct  correlation  function 
integrals,  but  significantly  affects  results  in  the  case  of  the  radial 
distribution  function  integrals.   The  previous  formalism  was  both  rigor- 
ous and  practical  for  properties  of  multicomponent  solutions  expressed 
in  terms  of  radial  distribution  functions  (assuming  the  constraint  to  be 
valid) ,  but  for  the  multisalt  case  yielded  intractable  equations  when 
expressed  in  terms  of  direct  correlation  functions.   This  is  primarily 
because  the  multicomponent  charge  neutrality  conditions  were  mishandled. 

The  method  used  here  parallels  that  of  the  previous  chapter  and 
yields  equivalent  results  to  those  of  O'Connell  and  DeGance  for  the  sin- 
gle salt  case,  and  quite  simple  results  for  multisalt  systems.   In  addi- 
tion, this  method  does  not  require  that  the  constraints  on  the  radial 
distribution  function  integrals  be  valid,  although  this  situation  is 

addressed  further  in  Appendix  G. 
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Material  Balance  Relationships  among  Salts,  Solvents  and  Ions 
In  the  case  of  strong  (completely  dissociating)  electrolytes  a  more 
compact,  but  equivalent,  set  of  matrices  are  used  in  the  material  balan- 
ces and  other  relations  than  were  used  in  the  general  formalism.   This 
set  is  used  because  of  the  sparseness  of  the  W  matrix  and  the  fact  that, 
with  the  salts  as  reference  components  (which  is  so  in  this  chapter) , 
the  K  matrix  is  a  zero  matrix.   The  new  matrices  are  defined  below. 

We  look  at  a  system  of  strong  electrolytes  dissolved  in  n  solvents, 
where  the  total  number  of  salts  plus  solvents  is  n   and  the  total  number 

of  ions  plus  solvents  is  n  (n  >  n   if  independent  salts  are  used  -  see 
Appendix  H) .   The  n  components  will  be  denoted  by  greek  subscripts  (a, 
B,...)  and  the  n  species  by  i,j,k,...  . 

We  construct  a  matrix  V^  defined: 
for  solvents 


V.   =  6. 
la    la 


i  =  1,.  .  .  ,s 
a  =  1 , . . .  ,s 


(4-la) 


and  for  electrolytes 


the  number  of  ions 

V.      =   of   type   i   in   a  molecule 
ict 

of  salt  a. 


i  =  s+1 , . . . ,n 
^   =  s+1 , . . . ,n 


(A-lb) 


with  all  other  elements  being  zeroes.   In  the  notation  of  Chapter  3,  the 
V  matrix  can  be  written 


V  = 


I  0 
0  V 

0  V, 


D 


(4-2) 


Such  a  system  has  n-n   constraints  on  its  composition  in  terms  of  ions 
and  solvents,  and  they  may  be  represented  by  an  n  by  n-n  matrix  z^ 
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having  the  property 


which  gives 


T 
V  z  =  0 


(4-3a) 


0 


(4-3b) 


The  elements  of  z^  include  the  ionic  charges  as  well  as  any  additional 
stoichiometric  constraints.   As  these  elements  are  not  functions  of  com- 
position, the  z  matrix  corresponds  to  the  Z  matrix  defined  previously. 
~  — o  '^  ■' 

In  order  to  perform  the  desired  projections,  two  more  nonsquare 
matrices,  u  and  j_,    are  defined  such  that 


V  u  =  J  , 


T 
u  z  =  0 


(4-4) 


^T  =  L  ^ 


TIl   ?^  0 


(4-5) 


where  J  is  an  n  xn   identity  matrix;   (J).. 
-        o   o  -  ij 


6 . .  ,  i, j  =  1 , . . . ,n  . 


The  _u  matrix  spans  the  same  nQ-dimensional  subspace  of  R  as  V  and  is 
orthogonal  to  z^  also.   The  form  of  _u  will  be  derived  below.   The  ^ 

matrix  is  a  very  simple  matrix,  but  is  not  orthogonal  to  z^.   It  has 
the  form 


I     0 


(4-6) 


and  from  equations  (4-5)  and  (4-6)  we  know  that 


u  =  y^  +  zb 


(4-7) 
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where  b  is  an  (n-n  )Xn  coefficient  matrix,  the  details  of  which  are 


not  important. 

Finally,  we  define  the  nXn  matrix 


T   =  vu^  ^  uv^   =  P^  (4-8) 


where  P^  is  the  orthogonal  projector  for  the  component  subspace,  projec- 
ting perpendicular  to  the  kernel  and  parallel  to  the  image  of  \;.   A 
property  of  P^  is  that  for  any  matrix  M  with  a  kernel  that  contains  all 
the  vectors  of  z_,  i.e.  Mz^  =  0^,  then 

MP  =  M  (4-9) 

Equations  (4-8)  and  (4-9)  also  imply  that 

P^  =  P  (4-lOa) 

Pu  =  u  (4-lOb) 

PV  =  V  (4-lOc) 

Pz  =  0  (4-lOd) 

and  with  equation  (4-7) 

PX  =  Z  (^-11) 

All  of  these  properties  of  P^  will  be  used  below. 

The  necessary  matrices  are  easily  constructed  from  the  known  \^ 
matrix.   Since  u_   and  V_   span  the  same  subspace  of  the  n-space,  each 

column  of  ii  must  be  a  linear  combination  of  the  columns  of  V.   In  all, 

2 
the  construction  of  u  from  V  requires  n   coefficients  which  make  up 

—      —  o 

the  n  xn  matrix,  c,  such  that 
o   o         — 

u  =  vc  (4-12) 
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Using  this  identification  in  equation  (A-4)  gives 


T     T 

V  u  =  V  VC  =  J 


(4-13a) 


or 


C  =  (v  v) 


(4-13b) 


This   identification  leads  to 


T   -1 
u  =  v(v  v) 


(4-14a) 


or 


u  = 


0    Vj^(v^D  +  V^Vj)-^ 

/  T    ,   T   ,-1 
0    Vj(Vj^D  +  VjVj) 


(4-14b) 


and 


/  T  ,-1  T 
P  =  v(v  V)   V 


(4-15a) 


or 


P  = 


1 
0 

0 


T     T   -1  T 


,  T   ,  T   ,-1  T 

V  (V  V  +V  V  )   V 


T     T    -1  T 


,  T   ,  T   .-IT 


(4-15b) 


The  material  balance  relationship  between  the  species  and  compo- 
nents is 


N  =  VN 


(4-16a) 


where  N   is  the  vector  of  moles  of  solvents  and  salts  (components).   The 


T   —  1       —  1 
v(v  v)   =  (v)    iff  V.  is  square  (not  possible  here) 
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corresponding  relations  from  species  to  components  are 


T      T 
N   =  u  N  =  y  N 
— o '^ 


(A-16b) 


As  before,  the  matrix  expression  relating  species  composition 
fluctuations  and  radial  distribution  function  Integrals  Is 


A  =  X  +  XHX 


(4-17) 


where  the  matrices  A,X  and  tl  have  been  defined  previously.   The  species 
composition  fluctuation  matrix  restricted  to  n   degrees  of  freedom,  A_, 
is  given  by 


T 

T    —I 
^  a 


9y 


oa 


T,V,y 


oY?^a 


A(X), 


where  the  elements  of  A   are  given  by 


(4-18a) 


^^)ae 


kT       oa 


3m 


o3 


9N 


kT 

N 


o3 


9V1 


oa 


T,V,y 


oy^a 


T,V,p 


M) 


6a 


oyt^B 


a, 3  =   1, . .  .  ,n 


(4-18b) 


Equations  (4-18)  can  be  derived  in  a  manner  analogous  to  equation  (3-37) 
using  the  methods  of  that  chapter.   Here,  however,  the  matrices  are 
slightly  different  (as  shown  above)  and  we  explicitly  define  the  salt 
chemical  potentials  by 


^   =  V  y  (4-19) 

where  JJ  is  the  vector  of  chemical  potentials  of  solvents  and  ions.   It 
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should  be  noted  that  although  this  is  somewhat  different  from  equations 
(2-4)  and  (2-5),  the  procedure  followed  in  Chapter  3  still  applies. 

As  in  O'Connell  and  DeGance  (1975),  A_   is  nonsingular  and  can  be 
inverted  to  give 


M-') 


N   oa 


a6   kT  8N 


o6 


N  ^^oB 

kT  9N 

oa 


T,V,N 


0Y7^3 


(A"') 6a 


(4-20) 


T,V,N 


oy^a 


The  radial  distribution  function  integrals  are  related  to  A_   by 

T 


T  .  ^i^I 


X  + 


''K 


T,V    J 


(X  +  XHX)y 


(4-21) 


This  equation  is  much  less  formidable  than  it  seems,  since 


T,V   J 


(X  +  XHX)z  =  0 


(E-11) 


or 


9y 


''^ 


=  -x'^(X  +  XHX)z(z'^(X  +  XHX)z}  ^ 


T,V 


(4-22) 


This  allows  us  to  rearrange  equation  (4-21)  to 


1 


A  =  T^l.   -  (X  +  XHX)z[z^(X  +  XHX)z]  ^z^}(X  +  XHX)y      (4-23) 
Of  course,  A        is  just  the  inverse  of  the  right  hand  side  of  (4-23). 


See  Appendix  G  for  certain  simplifications, 
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The  set  of  equations  that  equation  (4-23)  represents  can  be  solved 
by  simple  quadrature,  as  opposed  to  the  more  complex  method  required  to 
solve  the  system  of  partial  differential  equations  given  in  equation 
(3-57).   However,  even  these  quadratures  are  involved  due  to  the  mul- 
tiple matrix  operations  (especially  inversions)  which  need  performed  at 
every  quadrature  point.   This  difficulty  can  be  alleviated  if  the  dimen- 
sionality of  the  matrices  are  small  enough  to  allow  doing  the  matrix 
manipulations  analytically.   An  example  of  the  analytical  method  is 
given  below.   Thus,  as  before,  the  direct  correlation  approach  should 
be  of  greater  practical  value.   It  has  been  shown  (Gubbins  and  O'Connell, 
1974;  Brelvi  and  0' Connell , 1973, 1975)  that  direct  correlation  function 
integrals  are  insensitive  to  angle-dependent  forces  and  yield  excellent 
correlation  and  prediction  of  the  properties  of  gases  dissolved  in 
various  liquids,  including  water  (Mathias  and  O'Connell,  1980a, b). 

Formulation  in  Terms  of  Direct  Correlation  Functions 
O'Connell  and  DeGance  (1975)  showed  how  expressions  for  single 
salts  which  were  simple  functions  of  nondivergent  direct  correlation 
function  integrals  could  be  attained.   However,  due  to  incomplete 
expressions  for  charge  neutrality  and  composition  constraints,  this 
simplicity  disappeared  in  the  multisalt  systems.   An  ancillary  effect 
of  this  oversight  was  that  some  terms  were  stated  to  be  included,  but, 
in  fact,  are  not.    The  following  results  are  considerably  simpler. 


On  pages  767  and  777  of  O'Connell  and  DeGance  (1975),  it  is  noted 
that  all  terms  of  order  (M-1)  and  less  in  the  divergent  integrals  are 
preserved.   In  fact,  only  the  first  order  terms  remain. 
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Recognizing  that  the  direct  correlation  function  integrals  of  ions 
are  divergent,  we  return  to  the  definitions  of  the  direct  correlation 
functions  themselves.   For  species  of  type  i  and  j 


(r^^)    ^h   (r^^)  -I 


ij 


k=l 


\k(^23^Pk^kj^'^13^'^^3    (^-2^) 


In  order  to  perforin  the  transformations  necessary  to  attain  the  desired 
matrix  relations,  both  c..(r   )  and  h..(r   )  must  be  integrable.   How- 
ever, c^.(r   )  has  a  long  range  limit  (O'Connell  and  DeGance,  1975)  of 

2 

which  has  a  divergent  integral  over  volume.   Here,  e  is  the  charge  on 
an  electron  and  e,the  bulk  dielectric  constant  of  the  solvent  mixture. 
The  presence  of  the  z.z.  factor  in  equation  (A-25)  causes  these  portions 
of  the  direct  correlation  functions  to  lie  in  the  kernel  of  the  projec- 
tion operator,  P^.   Therefore,  these  divergences  can  be  removed  in  the 
process  of  projecting  the  species  composition  fluctuation  matrix.   We 
now  write  c..(r)  as  the  sura  of  a  nondivergent  part,  c..(r),  and  the 
divergent,  long-range  limit,  giving 


^ij(^)  =^ij(^>^Sf 


(4-26) 


Substituting  equation  (4-26)  into  (4-24)  and  operating  from  the  right 
with  P^  gives 


n   z .  e 


JAi^'n^^^h^  ^  J  ,  Tk^  ^j(^>jk  =  ^  ,  "ij(^i2)^^)ik 
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1=1 


hi£(^23^P?, 


n   z^e 


.1,  4j(^13)(^V-.J^^^,(Z)j. 


dr. 


(4-27) 
Since  the  vector  of  species  charges  (the  z.'s)  is  in  the  space  spanned 
by  z^,  we  have  by  equation  (A-lOd) 


n 


(4-28) 


Use  of  this  relation  in  equation  (4-27)  simplifies  it  tc 


n   n 


,!,^°,<^Iz)®,.  =  ^I/.j<^nH2V-j,j^ 


^£^^23^  P£ 


X  c^.(r^3)dr3 


(P) 


jk 


(4-29a) 


or,  after  Fourier  transforming, 


C°P  =  HP  -  HXC°P 


(4-29b) 


where 


(C°) 


IJ 


N 
V 


c^j(r;p,T)  dr 


(4-29c) 


Pre-multiplying  equation  (4-29b)  by  X  and  rearranging  gives,  with 
equation  (4-17) 

XHP  =  AC°P 


-u 


Adding  P  to  both  sides  and  using  AX   P.  =  £  +  XHP,  we  find 


A(X  ^  -  C°)P  =  P 


(4-30) 


which  is  nearly  identical  to  equation  (3-52) 
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We  project  to  component  variables  by  first  pre-multiplying  equation 


(4-30)  by 


T,V 


,  the  right  hand  side  simplifying  via 


equations  (4-lOd)  and  (4-11). 


T,V 


A(X  ^  -  C°)P  =  u^ 


(4-31) 


Equation  (E-11)  makes  it  clear  that 


T 

T  ,  ^% 


X  + 


^Ho 


T,V    J 


A   has  the 


required  characteristics  of  matrix  M  in  equation  (4-9).   Thus  we  know 


T 


T,V 


AP 


T 


T,V 


(4-32) 


Substituting  this  into  equation  (4-31)  and  using  (4-8)  gives 


T 


T,V 


T   -1 

AuV  (X 


C°)P  =  u^ 


(4-33) 


Post-multiplying  by  V^  and  successively  using  equations  (4-lOa) , (4-11) 
and  (4-18),  equation  (4-33)  becomes 


Av^iY.   ^   -   C")v  =  J 


(4-34a) 


or 


A-'   =  v^(x-^  -  cO)v 


(4-34h) 
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Equation  (4-34b)  provides  a  very  straightforward  relationship 
between  the  desired  macroscopic  salt  properties  (elements  of  ^  )  and 
the  molecular  direct  correlation  function  integrals  for  the  ions  and 
solvents.   Note  that  both  V  and  X   are  both  simple  and  well-known 
matrices  and  that  all  solvent-solvent  relationships  are  left  totally- 
unaffected  by  these  transformations.   This  result  is  identical  to 
equation  (3-77),  as  required. 


Thermodynamic  Properties 
The  general  matrix  relationships  between  measurable   thermodynamic 
properties  and  the  correlation  function  integrals  are  summarized  here. 
We  choose  to  work  with  the  mean  ionic  salt  properties,  such  as  the  mean 
ionic  activity  coefficient  on  the  mole  fraction  scale. 


£n  Y^   =  — 

±a   V 
a 


y  V.  £n  y. 
.  ^  la  'i 
i=s+l 


V  kT 
a 


y   V.  (y.-y.  -  kT  £n  N./N) 
i=s+l 


V  £n  Y^ 
a    ±a 


y  -  y 
o    og 

kT 


y    V.   £n 
i=s+l 


y  V.  N  /N 

1  ly  oy 


(4-35) 


with 


n 

V  =  y  V. 

a   .    ,      la 
1=1 


a=s+l , . . . ,n 


Here  we  have  tacitly  assumed  that  y   =  vy  ,  which  are  the  usual  standard 


states  chosen  for  salt  solutions.   Differentiating  with  respect  to  N 
yields  (Friedman  and  Ramanathan,  1970) 


oS 


9£nY 


Nv 


±a 


a  9N 


o3 


T,P,N 


oy/6 


N  ^%a 

kT  9N  o 
o3 


T,P,N 


n  V.  v.„ 

I     ^«^  +  ^  ^ 

i=s+l   ^      "  ^ 


(4-36) 


oy/6 
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where, 


X.  =  N./N 

X      1 


If  equation  (4-36)  is  restricted  to  salts  of  two  different  ions  (+  and 
-)  and  to  systems  without  ions  common  to  different  salts,  it  becomes 
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The  partial  molar  volumes  are,  for  solvents 
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and  for  salts 
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The  partial  molar  volumes  are  related  to  the  A_   matrix  by  the 
relation 
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where  o        =  N   /V,  k„  is  the  isothermal  compressibility  and  i  is  a 
oa    oa     T  — 

column  vector  of  n  ones.   Also, 


3P/RT 


3p 


1        T   -1  T 
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(4-AO) 


where  p  Is  the  molar  density  for  the  ions  plus  solvents,  p  =  N/V. 
Further, 

G  =  A~^  -  V  -   v'^X~^v  +  v^Iv  (4-41) 


where 
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(4-43) 


and  2.  is  an  nxn  matrix  of  ones. 

All  of  the  desired  properties  can  then  be  obtained  from  correla- 
tion function  integrals  using  equations  (4-34)  and  (4-39)  to  (4-41). 
The  expressions  in  terms  of  the  (H) . .  are  complex  (O'Connell  and  De- 
Gance,1975;  Friedman  and  Ramanathan,  1970)  whereas  they  are  simple  when 
the  ( C^  )  .  .  are  used. 
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-  [i-(c°),,][i-(c°).^]} 


(4-46) 


These  results  are  in  the  same  form  as  for  nonelectrolytes  (Mathias  and 
O'Connell,  1980a, b).   Equations  (4-45)  and  (4-46)  are  much  less  formi- 
dable than  they  appear  due  to  the  sparseness  of  \)  and  cancellations. 
For  salts  of  two  different  ions,  the  double  summation  over  i  and  j  leads 
to,  at  most,  four  terras. 

Since  they  can  be  modeled  well  in  nonelectrolyte  solutions,  we  also 
believe  it  is  useful  to  work  with  fluctuation  matrices  in  terms  of  den- 
sity rather  than  moles.   We  know  that  for  both  salts  and  solvents. 
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Defining  an  activity  coefficient,  y,  based  on  concentration  or  molar 
density. 
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We  can  write 
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and 
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Note  that  for  a=  salt,  we  write  V  £ny    instead  of  £ny   .   Previous 


a   ±a 


oa 


success  (Mathias  and  O'Connell,  1980a, b)  in  modeling  (C   ) . .    for  non- 
electrolytes  may  indicate  that  the  same  may  be  done  for  electrolyte 
systems.   In  this  case,  the  expressions  should  include  short-range, 
infinite  dilution  and  Debye-Hueckel  contributions.   However,  these 
are  ion-specific,  not  salt-specific. 

Single-Salt,  Single-Solvent  System 
The  simplest  example  of  salt   solutions  yields  results  obtained 
by  less  elaborate  methods.   However,  for  illustrative  purposes,  we  will 
use  the  above  method  to  derive  these  expressions.   In  this  system,  the 
solvent  is  indexed  1  and  the  salt,  with  ions  +  and  -,  is  indexed  2. 
Since  the  choice  is  arbitrary,  we  chose  +  as  the  dependent  and  -  as  the 
independent  ions.   This  choice  gives,  for  equations  (4-2)  and  (4-6) 
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The  z^   matrix  is,  from  equation  (4-3b)  , 
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If  we  multiply  z^  by  z  ,  it  becomes 
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which  shows  that  z_   is  indeed  parallel  to  the  vector  of  ionic  charges 
as  required  in  equation  (A-27) . 

First,  we  give  matrix^  in  terms  of  (H) . . .   Equation  (4-23)  gives 
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(4-53a) 
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and 


V  =  V  +  v_ 


(4-54) 


Interchanging  the  +  and  -  subscripts  does  not  affect  the  results.   This 
must  be  true  since  the  choice  of  dependent   and  independent  ion  was 
arbitrary.   The  A_       elements  can  be  determined  by  inversion. 

For  the  direct  correlation  function  integrals,  the  elements  of 
A_       can  be  determined  directly  from  equation  (4-34b).   From  equation 
(4-20) ,  we  now  find 
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v^2/^^RT=v{x^^(l-C°2)  +Vx^2^1-C°2)} 


(4-60) 
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The  direct  correlation  function  integrals  can  be  related  to  infin- 
ite dilution  values  and  the  Debye-Hueckel  limiting  law  expressions. 
The  correct  expressions  are 
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with  e,  the  solvent  dielectric  constant,  N  ,  Avagadro's  number  and  e, 
1  A 

the  unit  charge.   Then  (Redlich  and  Meyer,  1964) 


9£nY^ 


v  „  -  v  „  =  VRT  ^^ 
o2    o2        dP 


T  N    N 
'  ol'  o2 


(4-65a) 


1  S  (p  ,x  -)^  +  0(x  J 

2  V   ol  o2        ol 


(4-65b) 


where 


'9£ne, 


S  = 


S  RT 

Y 


9P 


(4-66) 


80 


and  K.  is  the  solvent  isothermal  compressibility.  The  solution  isother- 
mal compressibility  defines  the  apparent  molar  compressibility,  4),  ,(Guc- 
ker,  1933). 
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(General  expressions  for  the  multisalt  case  are  found  in  Appendix  I.) 
These  equations  yield 
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Double-Salt,  Single-Solvent  Case 
Assuming  the  two  salts  contain  different  ions,  we  write  the  v^,  ^ 
and  z^  matrices,  which  are  simple  extensions  of  the  one-salt  case.   Again 
we  choose  the  positive  ions  as  dependent,  and  the  negative  ions  as  in- 
dependent. 
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we 


Multiplying  the  first  column  of  z^  by  z   ,  and  the  second  by  z^_, 
again  get  the  matrix  of  ionic  charges  for  salts  2  and  3.   The  A_  matrix, 
in  terms  of   total  correlation  function  integrals,  becomes  quite  com- 
plex.  In  terms  of  component  matrices,  it  is 
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(4-73) 

The  elements  of  A        in  terms  of  the  total  correlation  function 
integrals  must  be  determined  by  carrying  out  the  matrix  operations  in 
equation  (4-73)  and  then  inverting  the  result.   It  is  clear  that  this 
would  be  a  very  complex  expression. 

However,  equations  (4-34b)  and  (4-72)  allow  us  to  write  the  ele- 
ments  of  A        easily,  in  terms  of  the  direct  correlation  function  integ- 
rals.  These  elements  are: 
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where  the  C  r,    (ct»3=l,2)  above  are  defined  in  equations  (4-56)  and  (4-57) 
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V2V3 


(4-78) 

In  the  above  equations,  the  other  derivatives  of  the  chemical  potential 
may  be  found  by  simply  interchanging  subscripts  2  and  3  and  defining 
C°   and  C°   in  a  manner  analogous  to  C   and  C„„.   Expressions  for  the 
thermodynamic  properties  are  the  same  as  those  of  Friedman  and  Ramana- 
than  only  when  the  constraint 


(X  +  XHX)z  =  0 


(4-79) 


is  assumed  (for  the  results  in  this  case,  see  Appendix  G) .   The  expres- 
sions in  terms  of  the  total  correlations  quickly  grow  cumbersome,  but 
the  expressions  in  terms  of  direct  correlations  remain  unchanged.   Since 


the  matrix  in  terms  of  the  (H)  . ..  must  be  inverted,  and  because  of  the 
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complexity  of  the  relationship  between  A_   and  the  (_H)  .  .  ,  we  prefer  to 

focus  on  the  expressions  using  direct  correlation  function  integrals, 

which  avoid  both  complexities.  Using  equations  (A-44)  to  (4-46),  we 
find 
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Expressions  involving  salt  3  can  be  obtained  from  those  in  equations 
(4-81)  to  (4-83)  by  interchanging  subscripts  2  and  3. 

The  preceding  examples  show  that  it  is  advantageous  to  write  the 
thermodynamic  properties  in  terms  of  salt-salt  and  salt-solvent  func- 
tions, taking  linear  combinations  of  the  ionic  direct  correlation  func- 
tion integrals.   If  we  note  that 
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equations  (4-44)  through  (4-46)  can  be  transformed  into  these  salt 
functions  in  the  completely  general  case.   We  have 
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and  V   =  1  if  a  is  a  solvent.   Similarly,  equation  (4-50a)  becomes 
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Single-Solvent,  Common  Ion  Case;  Harned's  Rule 
The  equations  change  somewhat  when  a  two-salt  system  has  a  common 
ion,  e.g.,  the  anion,  now  denoted  -.   Then 


C°2  =  [V^^(C°)2^2++  2V2^V2_(C°)2^_+  v^_(C°)__]/v2 


(4-92) 


C°3  =  [V2^V3^(C°)2^3^+  V2^V3_(C°)2^+  V2_V3^(C°)  3^_+  V2_V3_(C°)_J 


^  V2V3 


(4-93) 


C33  =  [V3+(C°)3_^3^+  2V3^V3_(C°)3^_+  v2_(c°)_J/v^  (4-94) 

and  the  number  of  different  direct  correlation  functions  has  been  re- 
duced from  fifteen  to  ten.   We  choose  the  two  positive  ions  as  the 
dependent  species. 

The  \;,  Y.   ^ri'i  ^  matrices  for  this  sytem  are 
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and 
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■V/^3+ 


(4-96) 
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where  again,  multiplication  by  z_  gives  the  vector  of  ionic  charges. 
Since  the  C  ^'s  use  the  V  matrix  directly,  whereas  the  (H)..'s 
must  be  transformed  by  z_,    inverted,  further  transformed  by  combinations 
of  X,  H,  z^,  and  _y^,  then  inverted  again,  it  is  obvious  that  thermody- 
namic properties  in  terms  of  the  direct  correlation  functions  are 
much  simpler  mathematically. 

Harned's  rule  for  the  mean  ionic  activity  coefficients  states 

that  the  logarithm  of  the  mean  ionic  molal  activity  coefficient,  y   , 

m 

of  one  salt  varies  proportionally  to  the  molality  of  the  other  salt 
when  the  total  ionic  strength,  I,  is  kept  constant.   That  is 
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where  v   -.  and  I  are  defined  by 
m±z 

I  =  (^2+^2+  "^  ^2-^-^"'2  "^  ^^3+4+  "^  ^3-^?^'"3   ,       (4-98b) 

m   is  the  molality  of  salt  a  and  m,„  is  the  mean  ionic  molality  of 
a  ±2 


salt  2. 


V„   V   1/v 
m^2  =  (n'2+  '"_   )  ^^-'^"^^ 


From  the  definitions  of  density-  and  molality-based  activity 
coefficients  in  equations  (4-48)  and  (4-98a) ,  equation  (4-97)  can  be 
rewritten 
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Next,  we  use  equation  (4-50)  to  write 
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(4-101) 

The  formulae  for  the  necessary  direct  correlation  function  integ- 
rals can  be  obtained  from  Appendix  I.   They  are 
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(4-103) 


2  ^     2 
where   q   =  V  , z  ,  +  v   z 

a    a+  a+    a-  a- 
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The  expressions  for  1-C„„  can  be  obtained  by  replacing  subscript 

3  by  2  in  equation  (4-103),  except  in  the  term  (x  „q„  +  x  ^q^) ,  where 

ol   I  o3  3 

the  3  remains.   The  paramter  B   will  contain  short  range  contributions 
from  ions  in  salts  a  and  3  and  long  range  terms  in  a  universal  form 
from  the  salts. 

Equations  (4-102)  and  (4-103)  are  substituted  into  (4-101)  and 
considerable  cancellation  of  terms  occurs  (including  all  terms  propor- 
tional to  the  inverse  of  composition)  leaving  a  series  expression  for 
the  Harned  coefficients,  in  terms  of  molality. 
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(4-104) 

where  M  is  the  molecular  weight  of  the  solvent  (component  1) .   It  is 
of  interest  to  note  that  the  Harned  coefficient  is  a  weak  function  of 
total  ionic  strength,  which  is  constant  for  a  given  system,  as  well  as 
of  m  ,  which  varies. 

For  a  1:1  electrolyte,  the  zeroth-order  approximation  to  a   is 


°'23   1000  ^^22   ^23 


lih'i^^li-^,^^^'^'"'^ 


(4-105) 


The  values  of  B  „  and  B   will  have  separate  contributions  from  the 
ions  and  will  depend  only  on  temperature.   Thus,  in  principle,  a 
should  be  correlatable  using  only  pure  salt,  pure  solvent  and  Infinite 


90 


dilution  data.   Conversely,  the  determination  of  Harned  coefficients 
will  aid  in  establishing  expressions  for  the  (C°) 


ij 


CHAPTER  5 
ELEMENTS  OF  A  MODEL  FOR  ELECTROLYTIC  SOLUTIONS 


Introduction 
The  previous  chapters  give  a  formal  basis  for  calculating  system 
properties  from  direct  correlation  function  integrals,  but  deal  very 
little  with  precisely  how  these  calculations  can  be  performed  or  how 
these  properties  can  be  correlated  or  predicted  using  direct  correla- 
tion functions.   Previous  work  (Mathias  and  O'Connell,  1979,  1980a, b) 
has  shown  that  properties  for  systems  of  non-"reactive"  supercritical 
components  dissolved  in  various  solvents  can  be  correlated  and  pre- 
dicted quite  well.   In  these  systems,  two  pure-component  parameters 
and  one  condition-independent  binary  parameter  were  sufficient  for  the 
accuracy  desired.   A  more  general  result  of  this  work  is  to  indicate 
that  the  use  of  direct  correlation  function  solution  theory  to  model 
solution  properties  is  most  appropriate  in  systems  where  the  compo- 
nents are  very  different  from  each  other  -  hence  the  very  good  results 
for  a  mixture  of  supercritical  and  a  subcritical   component.   This 
result  stems  from  the  mathematical  form  of  the  equations  which  amplify 
the  effect  of  modeling  errors  when  the  values  of  the  direct  correlation 
function  integrals  are  similar,  and  not  from  restrictions  in  the  theory 
itself.   With  this  knowledge  and  the  work  of  the  previous  chapters  In 
mind,  solutions  of  electrolytes  seem  to  be  a  natural  application  of 
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this  theory,  due  to  the  "reactive"  nature  of  the  electrolytic  dissocia- 
tions and  the  differences  in  the  physical  nature  of  the  solvents  and 
ions. 

This  chapter  concerns  the  basis  and  general  formulation  of  a 
simple  model  for  electrolytes,  and  how  this  model  can  be  used,  e.g.  in 
determining  activities  of  electrolytes  or  the  effect  of  salts  on  vapor 
liquid  equilibria.   The  model  is  not  complete,  but  a  few  preliminary 
results  are  given. 

The  Primitive  Ion  Model 

The  primitive  ion  model  approximates  an  ion  in  solution  as  a 
hard  sphere  with  an  electric  charge  at  its  center.   Thus,  an  ion  i  is 
completely  described  by  charge  parameter  z.  (number  of  elementary 
charges)  and  a  size  parameter  CT_  (hard  sphere  diameter)  .   There  are  a 
number  of  theories  based  on  this  view  of  an  ion.   The  simplest  is  the 
Debye-Hueckel  Theory,  which  looks  at  one  ion  in  a  uniform  field  gen- 
erated by  all  the  remaining  ions,  and  in  its  original  form  considers 
the  ion's  diameter  to  be  zero.   This  theory  gives  exact  results  for 
very  dilute  solutions  (e.g.  10   molar  aqueous  solutions  for  ions  with 
single  charges).   A  more  recent  approximation  to  a  solution  of  primi- 
tive ions  is  known  as  the  Mean  Spherical  Model  (Lebowitz  and  Percus, 
1966).   It  has  been  solved  for  both  the  "restricted"  case,  where  all 
ions  have  equal  sizes  and  charges  (Waisman  and  Lebowitz,  1970),  and  in 
the  general  case  of  arbitrary  sizes  and  charges  (Blum,  1975) .   Un- 
fortunately, this  model  requires  tedious  iterative   calculations  in 
all  but  the  restricted  case,  and  does  not  give  the  correct  results  in 
the  limit  as  the  salt  concentration  goes  to  zero  (i.e.  Debye-Hueckel 
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results) .   This  model  gives  good  results  for  solutions  of  moderate 
concentrations  in  both  the  general  and  restricted  cases  (Grigera  and 
Blum,  1976;  Triolo  et  al. ,  1976). 

The  model  proposed  here,  also  using  the  primitive  ion  basis,  is 
easier  to  use  than  the  Mean  Spherical  Model,  but  should  give  reasonable 
accuracy  in  moderately  concentrated  solutions  of  electrolytes,  as  well 
as  go  to  the  Debye-Hueckel  limits  as  the  salt  concentration  vanishes. 
Work  by  Long  and  McDevit  (1952),  Shoor  and  Gubbins  (1969)  and  Lucas 
(1969)  indicate  that  the  major  portion  of  the  activities  of  mixtures  of 
solvents  and  electrolytes  stems  from  the  volumetric,  or  size,  effects 
of  the  individual  components,  rather  than  the  electrostatic  effects. 
For  this  reason,  we  choose  to  model  the  direct  correlation  function 
integral  as  the  sum  of  a  hard-sphere  term  plus  a  deviation,  or  pertur- 
bation, which  will  largely  be  from  the  Debye-Hueckel  theory.   This 
parallels  the  approach  of  Mathias  and  O'Connell  (1979),  where  the  per- 
turbation to  the  hard-sphere  integral  was  a  second  virial  coefficient 
term.   We  use 


C.  =  p 
11 


c  .  .  (£)  dx  =  p 


c.  .  (r)dr  +  p 
ij  ^  " 


[c..(r)  -  cl'!(r)]dr      (5-1) 


where  the  hard-sphere  term  can  be  obtained  from  differentiation  of  the 
mixture  form  of  the  Carnahan-Starling  equation  of  state  for  hard-spheres 
(Mansoori  et  al. ,  1971).   The  hard-sphere  portions  of  the  direct  corre- 
lation function  integrals  and  activity  coefficients  are  given  in  Ap- 
pendix I.   The  perturbation  term  comes  from  the  electrostatic  contribu- 
tions and  includes  the  Debye-Hueckel  limits  given  in  equations  (4-69) 
through  (A-70),  as  well  as  terms  from  the  Mean  Spherical  Model. 
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Calculation  of  Solution  Properties 
Before  the  model  for  these  direct  correlation  function  integrals 
is  examined  in  more  detail,  it  is  important  to  understand  the  process 
used  in  carrying  out  an  actual  calculation.   The  C.'s  are  written  as 
functions  of  T  and  p_  and,    as  shown  in  equations  (4-40)  through  (4-50), 
are  related  to  the  following  density  derivatives  of  thermodynamic  prop- 
erties : 

^'/-^        =   I   V.    I     ^[1-    (C°)..]  (5-2a) 

'  oYT^a 


3p 
oa 


and 


9£ny 


^PoB 


-  ^ap/Po  (^-2b) 


''W3 


In  a  typical  problem,  the  activity  coefficients  of  the  components 

of  a  solution  are  desired  at  a  particular  state  —  usually  defined  by 

the  variables  T  ,  P   and  x  .   Since  the  (C  )..'s  are  related  to  deriv- 

-o  -  ij 

atives  of  the  quantities  of  interest,  it  is  obvious  that  we  need  to 

know  the  complete  set  of  activity  coefficients  at  some  reference  state 

r   r      r 
given  by  T  ,  P   and  x    .   (Actually,  values  of  the  activity  coefficients 

at  reference  conditions  are  needed  only  for  those  components  whose  ac- 
tivity coefficients  are  to  be  calculated  at  T  ,  P   and  x  ) .   Further- 

— o 

more,  since  equatons  (5-2)  give  only  density  derivatives,  the  reference 

r   f 
state  temperature  must  equal  the  final  temperature  (T  =T  ) ,  and  the 

reference  state  density  must  be  known  (p  ) . 
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The  simplest  procedure  for  the  calculation  is  to  make  a  guess 

,    f      f      f      i 
of  the  density  at  the  final  state  (T  ,P  ,p  ,x   ) ,    noting  that  the  com- 
ponent densities,  the  p   's,  are  given  by  p   =  x   p.   Then  equation 
(5-2a)  is  integrated  successively  for  each  component  (i.e.  from  state 

{T  ,p    p   ,...}  to  state  {T^,p   p   ,..  }to  state  {T^,  p   .  p   ,...} 
oi   oz  qI      oz  ol   o2 

and  so  forth).   This  gives  a  pressure  P'  which  corresponds  to  the  de- 
sired mole  fractions,  but  probably  incorrect  total  density.   We  then 
use  equation  (4-44) ,  which  can  be  rewritten 
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Equation  (5-3)  is  integrated  to  give 

f 


p   _  p. 


RT 


^   ^  ^on^o 
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oa'^og 
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C„g(3,X^.T)dp 


(5-4) 


Integration  of  equation  (5-2a)  can  be  performed  analytically  or  numer- 
ically -  in  either  case  it  is  a  simple  quadrature.   The  integration  in 
equation  (5-4)  can  also  be  done  either  numerically  or  analytically.   If 
done  numerically,  simple  quadrature  is  performed  until  the  value  of  the 
right  hand  side  equals  the  known  value  of  the  left  hand  side,  the  end 
point  of  the  integration  interval  being  p  ,  the  desired  quantity.   If 
done  analytically,  a  root  finding  procedure  for  p  must  be  employed. 
Because  these  are  simple  integrations  rather  than  solutions  to  differ- 
ential equations,  they  can  be  done  quickly  and  efficiently  on  a  com- 
puter. 
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Once  p   is  known,  equation  (5-2b)  can  be  integrated  using  the 
same  procedure  as  followed  with  equation  (5-2a) ,  except  that  now  the 
final  limit  is  kno\^m  with  certainty  and  further  integration  at  con- 
stant composition  is  unnecessary.   These  integrations  can  also  be  per- 
formed either  numerically  or  analytically  (see  Appendix  J   for  inte- 
grated form  of  equation  5-2b) . 

The  flexibility  that  this  process  has  due  to  the  use  of  a  refer- 
ence state  should  not  be  overlooked.   For  a  system  of  salts  in  a 
single  solvent,  the  use  of  the  pure  solvent  as  the  reference  state, 
and  correspondingly,  considering  infinite  dilution  to  be  the  refer- 
ence state  of  the  salts,  only  the  solvent  density  at  an  arbitrary  ref- 
erence pressure  (the  temperature  must,  of  course,  be  the  same  as  the 
final  temperature)  is  needed  for  reference  data.   A  more  important 
example  of  judicial  selection  of  the  reference  state  occurs  in  vapor- 
liquid  equilibria  calculations.   Assuming  we  have  one  salt  in  some 
liquid  mixture,  the  salt-free  mixture  with  the  same  mole  ratios  of 
non-electrolyte  components  can  be  viewed  as  a  single  pseudo-component, 
and  this  solvent  used  as  the  reference  state.   From  experiments,  or 
theories  of  non-electrolyte  solutions,  the  values  of  p   and  the  non- 
electrolyte  activity  coefficients  can  be  obtained.   This  procedure 
allows  for  interactions  between  the  solvent  constituents  to  be  taken 
into  account  by  more  accurate  models.   Finall}',  the  integration  of 
equations  (5-2)  need  only  be  performed  over  the  salt  number  density 
and  solvent  "pseudo-component"  density,  rather  than  the  entire  set  of 
densities.   This  procedure  allows  the  salt-effects  on  the  solution 
properties  to  be  determined  using  our  model,  but  does  not  restrict  the 
solvent-solvent  interactions  to  the   primitive  ion  restrictions. 
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Thus,  the  heavy  dependence  on  a  reference  state  has  several  ef- 
fects, both  positive  and  negative.   It  requires  accurate  density  and 
activity  data  at  the  reference  state,  but  allows  us  to  make  us  of 
any  available  data  that  might  minimize  the  length  of  the  integration 
path.   And,  it  allows  us  to  be  concerned  only  with  the  effects  of  the 
electrolytes  and  total  density,  and  so  make  use  of  other  methods  which 
might  be  more  appropriate  for  the  solvents  of  interest  when  accounting 
for  solvent-solvent  interactions. 

Electrostatic  Contributions  to  the  Direct  Correlation 

Function  Integrals 

Since  the  largest  contribution  to  the  effect  of  ions  in  solution 
properties  is  the  effect  of  the  ion  size,  accounted  for  by  the  hard 
sphere  portion  of  the  direct  correlation  function  integral,  the  elec- 
trostatic contributions  can  be  included  in  an  approximate  way  without 
a  great  loss  in  accuracy.   Therefore,  the  superposition  of  the  Debye- 
Hueckel  limits  on  the  hard  sphere  relation  is  considered  adequate  for 
the  electrostatic  effects  on  solvent-solvent  and  a  valid  approximation 
for  ion-ion  interactions.   Mathematically  speaking,  the  Debye-Hueckel 
limits  can  be  looked  at  as  truncations   of  half-power  expansions  in 
salt  density,  with  the  largest  exponents  being  tliree-halves,  one-half 
and  negative  one-half  respectively   for  solvent-solvent,  solvent-ion 
and  ion- ion  interactions  (equations  I-^a  tlirough  I-Ab).   Since  the 
largest  effect  is  formed  in  ion-ion  interactions,  it  is  necessary  to 
extend  that  series  to  higher  powers  of  the  salt  density.   It  is  here 
that  some  results  from  the  Mean  Spherical  Model  are  used. 
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Beginning  with  the  general  solution  of  the  Mean  Spherical  Model 
(Blum,  1975)  and  approximately  expanding  in  half-powers  of  the  salt 
density  gives  indications  that  the  electrostatic  contribution  to  the 
direct  correlation  function  integrals  can  be  written 


elec       „   ,   °°'  k   ^ 

-  C..     =2:  r  I     D..,  (a,n)K  (P2L)  (5-5) 


with  the  following  definitions  of  the  terms  on  the  right  hand  side: 


z.   =   charge   (in  units  of  elementary  charge)  on  ion  i, 


2     4TTe^ 
a   =  — r-zr  and  has  units  of  length, 
e  kT  ^ 

o 


2      y  2 

K        =   a  I   p .z .    =   inverse  debye  length. 


O.      =      radius  of  ion  i. 


n,   =   a  weighted  average  of  ion  sizes,  having  units  of 
(length)*^  and  defined  by 
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It  is  important  to  note  that  if  we  factor  out  the  total  ionic  density 

(    I        p.)  from  both  numerator  and  denominator  of  the  Hu's,  they  are  seen 

ions 
to  be  functions  of  the  ratios  of  the  salt  mole  fractions  among  them- 
selves.   That  is,  if  the  salts  in  a  solution  are  held  in  constant 
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proportion  to  each  other,  the  Hu's  remain  unchanged  if  the  total  salt 
density  is  changed.  Thus  x]  is  not  a  function  of  the  total  salt  con- 
centration, only  the  relative  concentrations. 

Finally,  we  must  define  the  functions  D-.,  (a,jn).   A  general  form 
for  these  functions  has  not  been  determined,  but  the  dependence  on  o; 
and  j}  has  been  determined  for  k  values  ranging  from  zero  to  four. 
These  are 
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,3,3,  ,  ^  ^  ,2  ^     2 


c^{3/2ti^(o^  +  op   +  o.OjCa^  +  a^)  +  3/2113(0.  +  o^)  +  3/An^n3} 


+  d  {2aja^  +  3n2(aJ  +  op   +  3/4112}  +  e^{a.  +  ^^  +  3/2n^} 


In  the  above  equations,  the  a.'s,  b.'s,  c.'s,  d,  and  e,  are  pure 
numbers  arising  from  multinominal  and  binominal  expansion  coefficients, 
and  are  very  complex.   The  Mean  Spherical  Model  is  only  an  approxima- 
tion to  the  ions'  behavior,  and  in  fact  gives  an  incorrect  result  at 
infinite  dilution.   Furthermore,  this  truncated  expansion  is  only  an 
approximation  to  the  low  density  (<1  Normal)  limit  of  this  model. 
So  rather  than  evaluating  these  constants,  we  consider  them  to  be  uni- 
versal parameters  to  be  fitted  from  data.   Furthermore,  in  the  k=.2ero 
term  (which  is  proportional  to  salt  density  to  the  minus  one-half 
power) ,  a   is  replaced  by  the  correct  Debye-Hueckel  limit  of  equation 
(I-4c). 

The  Form  of  the  Complete  Model 
The  individual  contributions  to  the  model  for  the  direct  cor- 
relation function  integrals  have  been  defined  in  the  preceding  sections 
and  Appendices  I  and  J.   Now  these  parts  must  be  combined  to  form  the 
total  model.   We  consider  a  system  of  a  mixed  solvent  pseudo-component 
denoted  by  the  subscript  s,  and  n   salts  denoted  by  the  subscripts  a 

and  3,  where  a,B  range  from  1  to  n  .   The  general  equations  for  the 

_o  t        . 
C  „  s  are  then 
aB 
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c°  =  1  -  c^j-  (1  -  c^r  -  c^''+-^)  I V 
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aS  \)   %)„   ^  ^   la  iR      11 
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,  2  2°° 

+  a  p  y  y  v.^v.nz.z;.   y  D.  .,  I 

a  3  k=l   -^ 


k-1 


(5-7c) 


In  these  equations  K     and  e  are  the  isothermal  compressibility  and 
dielectric  constant  of  the  mixed  solvent  (salt  free) ,   The  other  var- 
iables are  defined  in  the  sections  dealing  with  the  specific  portions 
of  the  direct  correlation  function  integrals.   The  superscripts  hs, 
hso,  hs,  and  hs£  are  hand  sphere  contributions  and  are  given  in  Ap- 
pendix J. 


Preliminary  Results 
Although  fitting  of  the  parameters  in  the  D..,  's  is  beyond  the 
scope  of  this  work,  some  preliminary,  approximate  calculations  were 
made  for  the  sodium  chloride-water  system.   These  calculations  were 
carried  out  by  fitting  the  three  direct  correlation  function  integrals 
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(water-water,  water-salt,  and  salt-salt)  to  compression  data  from 
Gibson  and  Loeffler  (1948)  and  the  activity  data  reported  by  Robinson 
and  Stokes  (1965).   Then  with  the  dielectric  constant  and  its  deriva- 
tives taken  from  the  data  of  Owen  et  al.  (1961),  the  water-water  and 
water-salt  integrals  were  recalculated  according  to  equations  (5-7) . 
Using  the  known  final  density  [instead  of  using  equations  (5-2a)  and 
(5-4)]  and  infinite  dilution  partial  molar  volumes  of  the  ions,  and 
the  above  C^g  values  ,   the  thermodynamic  properties  of  the  solution 
were  calculated.   The  results  are  given  in  Table  1  along  with  the  o, 
values. 

This  procedure  gives  calculated  values  of  reduced  compressibility 
within  five  per  cent  error,  reduced  partial  molar  volumes  within  nine 
per  cent  (one  per  cent  for  NaCl)  and  derivatives  of  the   activity  coef- 
ficient within  twenty  per  cent  at  concentrations  as  high  as  five  to 
six  molar.   The  errors  are  systematic,  and  may  be  correctable  without 
too  much  difficulty.   For  example,  although  using  the  above  procedure 
to  integrate  (5-2a)  or  (5-4)  may  lead  to  inaccurate  densities;  it  may 
lead  to  more  accurate   actvities  than  is  implied  here  because  consis- 
tent application  of  the  model  will  have  been  made.   Further  refine- 
ments to  the  model,  as  well  as  correlation  of  the  higher  terms  in  the 
ion-ion  direct  correlation  function  should  lessen  these  errors. 

As  hoped,  this  modeling  procedure  is  a  viable  method  for  correla- 
ting and  predicting  properties  of  salt  solutions.   Further  refinements 
are  needed,  of  course,  yet  this  beginning  appears  to  be  quite  promising, 
Further  work  on  improving  the  accuracy  is  presently  being  undertaken. 
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CHAPTER  6 
AN  EXPERIMENT  FOR  VOLUMETRIC  PROPERTIES  OF  SALT  SOLUTIONS 

AT  HIGH  PRESSURES 


Introduction 

A  molecule  in  a  solution,  held  at  a  constant  temperature  and  com- 
position, would  notice  changes  in  the  pressure  of  the  solution,  not  so 
much  by  the  force  exerted  on  the  boundary  of  the  solution,  as  by  the 
nearness  of  other  molecules  in  the  solution.   This  concept  leads  to 
the  notion  that  when  looking  at  solutions  in  terms  of  their  inter- 
molecular  interactions,  the  density  becomes  the  fundamental  variable, 
and  the  pressure  a  measured,  or  perhaps,  derived  property.   It  is  the 
fact  that  the  density  is  a  more  direct  measurement  of  the  interactions 
between  molecules,  that  it  is  of  great  interest  to  look  at  the  varia- 
tion of  thermodynamic  properties  along  isochoric,  rather  than  isobaric, 
surfaces  in  P-V-T-x  space.   Further  motivation  is  given  by  calculations 
by  O'Connell  (1980)  from  experimental  data  by  Gibson  and  Loeffler 
(1948) .   There  it  is  shown  that  the  mole  fraction  derivative  of  the 
salt  activity  coefficient  of  sodium  chloride  in  water  is  a  much  smoother 
and  weaker  function  of  composition  along  isochores  than  along  isobars. 

The  purpose  of  high  pressure  volumetric  measurements  of  salt 
solutions  is  to  be  able  to  look  at  activities  and  activity  coefficients 
of  the  components  of  such  solutions  on  these  surfaces  of  constant  den- 
sity, the  elevated  pressures  being  necessary  because  of  the  relative 
insensitivity  of  density  to  pressure  as  compared  to  its  sensitivity  to 
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composition  and  temperature.   For  example,  water  at  1  atm  and  25°C  has 

-4   -1 
a  coefficient  of  thermal  expansivity  of  2.57  x  10    K   ,  and  a  coef- 

-5    -1 
ficient  of  isothermal  compressibility  of  4.52  x  10   bar   .   These 

values  correspond  to  a  pressure  increase  of  approximately  5.7  bar  for 
every  one  degree  increase  in  temperature  along  the  line  of  constant 
density  in  the  P-T  plane.   Because  of  similarly  strong  composition  ef- 
fects on  density  (especially  at  lower  salt  concentrations),  one  can  see 
that  a  large  pressure  range  is  required  if  a  reasonable  temperature  and 
composition  range  is  used. 

This  chapter  deals  with  the  use  of  the  P-V-T-x  data  to  calculate 
activities  and  activity  coefficients,  as  well  as  the  operation  and 
design  of  the  apparatus  itself.   Since  the  apparatus  is  that  of  Grindley 
and  Lind  (1971)  and  Grindley  (1971)  with  only  minor  changes,  the  dimen- 
sions, details  of  calibrations  and  equipment  specifications  are  left  to 
those  references.   Only  a  more  general  description,  including  altera- 
tions, is  given  here. 

Thermodynamic  Basis 
Given  that  data  from  this  experiment  is  available  over  a  spectrum 
of  P-V-T-x  values  for  the  solution  of  interest,  we  now  look  at  the  cal- 
culation of  the  activities  and  activity  coefficients  on  a  constant  den- 
sity surface.   The  variation  of  activity  with  temperature  and  composi- 
tion cannot  be  derived  from  P-V-T-x  data  alone,  but  the  variation  with 
respect  to  P  at  constant  T  and  X  is  given  by 


9£na 
oa 


9P 
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=  [v^^(T,P,x)  -  v^^(T°,P°.x°)]/RT  (6-1) 

T,x 


where  the  superscript  zero  refers  to  the  reference  state  used  for  the 

activity.   Equation  (6-1)  can  also  be  used  for  the  activity  coefficient 

by  replacing  a   with  Y 
^       ^  *  oa      'oa 

The  partial  molar  volume,  v   ,  can  be  calculated  by  fitting  the 

oa 

P-V-T-x  data  and  differentiating.   For  data  in  this  form,  the  equation 


v   =  v  -   y   X 
oa       „f2   oa 
B^a 


3v 


3x  Q 
o3 


(6-2) 
^'^'""oY^B.a 


is  the  most  direct.   Once  all  the  mole  fraction  derivatives  are  deter- 
mined, all  the  partial  molar  volumes  can  be  easily  calculated.   Simi- 
larly, the  partial  molar  volume  at  the  standard  state  is  determined 
using  equation  (6-2)  with  data  from  around  the  standard  state. 

Of  course,  the  P-V-T-2C  data  only  gives  the  changes  in  activity 
as  the  pressure  changes.   So,  we  need  to  know  the  activity  of  the  com- 
ponents in  a  solution  at  the  temperature  and  composition  of  interest 
and  some  other  pressure.   Also,  we  need  to  supplement  our  measurements 
with  published  values  of  activities  at  normal  pressures,  which  are 
readily  available  for  many  systems.   Having  the  activity  of  component 
a  at  T,  P' ,  2£j  where  P'  is  the  pressure  from  the  outside  reference,  we 

can  calculate  £na   at  any  pressure  using  the  relation 
oa      -^  ° 
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Jlna  ^(T,P,x)=  Jlna   (T,P',x)  +- 


J-    r  -  D 


P' 


^\a  -   ^oa^dP  (6-3) 


as  long  as  the  v   's  are  available  at  T  and  X  over  the  complete  pres- 
sure  range  from  P'  to  P.   If  the  standard  state  for  the  activity  is 
a  fixed  pressure  P   and  not  the  system  pressure  (which  varies  from 
P'  to  P) ,  we  may  rewrite  equation  (6-3)  as 


-o 

)lna   (T,P,x)  -  in^      (T,P',x)  -  -^^— (P  -  P')  +  — 
oa     -       oa'     '-     Rx  RT 


fP  _ 

V   dP 
,   oa 
P' 


Equations  (6-1)  through  (6-3)  enable  the  activities  to  be  cal- 
culated over  the  entire  range  of  P-V-T-x  data,  assuming  that  they  are 
available  at  some  pressure  for  the  same  T-x  range.   Since  the  volumes 
are  measured  throughout  the  range  of  P-T-x,  the  isochoric  surfaces  can 
be  constructed  for  the  activity  in  the  space  of  the  variables  T  and  X. 

Experimental  Design 
The  apparatus  used  consists  of  two  major  components.   One  is  the 
piezometer  apparatus,  including  a  pressure  tube,  piezometer  tube,  and 
volume  measuring  devices  contained  in  a  constant  temperature  bath.   The 
other  is  the  pressure  plant.  Including  pressure  generating  and  mea- 
suring devices.   The  main  portions  of  these  components  are  sketched 
in  Figures  6  and  7   respectively.   These  figures  are  described  in  the 
following  sections. 
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Piezometer  Design 

The  piezometer  tube  Itself  Is  a  non-magnetic  (austenitic)  stain- 
less steel  tube  which  holds  the  sample.   One  end  of  the  tube  is  sealed 
with  a  magnetic  (martensitic)  stainless  steel  plug.   The  other  end  has 
an  identical  plug,  but  it  is  removable.   This  allows  a  snugly  fitting 
martensitic  stainless  steel  float  to  be  placed  against  the  test  liquid, 
and  for  the  remainder  of  the  tube  to  be  filled  with  mercury.   VJhen  the 
plug  is  in  position,  a  small  hole  in  the  piezometer  tube  wall  allows 
pressure  to  be  transmitted  to  the  mercury  inside  the  tube,  and  hence 
to  the  sample  itself.   The  mercury  acts  as  a  seal  to  prevent  any  leak- 
age around  the  float.   The  piezometer  tube  fits  inside  a  mercury  reser- 
voir tube,  and  then  the  entire  assembly  slides  into  the  pressure  tube. 
Here  the  mercury  contacts  the  oil  which  is  the  pressure  transmitting 
fluid  in  the  pressure  plant.   This  configuration  is  shown  in  Figure  6. 

The  pressure  tube  is  fixed  to  a  bracket  which  also  holds  the 
mechanical  portion  of  the  volume  measuring  apparatus.   This  portion 
consists  of  a  linear  variable  differential  transformer  (LVDT)  which 
is  attached  to  a  motor-driven,  lathe  lead  screw.   The  LVDT  detects 
the  location  of  the  upper  plug  and  the  float  in  the  piezometer  tube. 
The  lathe  lead  screw  then  allows  the  position  of  the  LVDT  to  be  known 
very  accurately.   A  counter  attached  to  the  motor  drive  is  used  to 
count  rotations  of  the  lead  screw,  which  can  be  converted  into  the 
length  of  the  sample  liquid  columns  in  the  piezometer  tube. 

The  LVDT  is  simply  three  coils,  the  primary  coil  through  which 
an  AC  signal  is  sent,  and  two  identical  secondary  coils  which  are 
arranged  symmetrically  around  the  primary  coil.   The  secondary  coils 
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PIEZOMETER     APPARATUS 


Figure  0.   Piezometer  apparatus  for  volume  determinations- 
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are  wired  such  that  the  total  voltage  across  the  two  coils  is  zero  when 
the  martensitic  stainless  steel  plug  or  float  is  centered  between  the 
secondary  coils.  By  adjusting  the  lead  screw,  the  output  from  the  LVDT 
can  be  nulled,  which  corresponds  to  locating  the  desired  position.  The 
difference  between  the  location  of  the  plug  and  that  of  the  float  gives 
the  length  of  the  liquid  column. 

Pressure  Plant 

Pressure  is  generated  in  the  apparatus  by  operation  of  a  hydraulic 
hand  pump  connected  to  the  low  pressure  side  of  a  pressure  intensifier. 
The  intensifier ' s  high  pressure  side  is  connected  to  the  piezometer 
apparatus,  where  pressures  in  excess  of  8000  bar  can  be  attained. 
Another  hand  pump,  the  reversing  pump,  is  used  to  return  the  pressure 
transmitting  piston  of  the  intensifier  to  its  initial  position  after 
the  pressurization  is  completed.   A  valve  with  its  valve  seat  drilled 
out  is  used  for  fine  pressure  adjustments.   For  pressure  measurements, 
the  pressure  plant  includes  a  300  bar  Heise  gauge  and  40,000  psia 
Heise  gauge  which  are  used  for  low  and  intermediate  pressures,  respec- 
tively.  Intermediate  and  high  pressures  are  determined  using  a  man- 
ganin  gauge  which  is  described  in  the  following  paragraph.   A  schematic 
of  the  entire  pressure  plant  is  given  in  Figure  7. 

The  manganin  gauge  measures  pressure  by  measuring  the  change  in 
resistance  of  a  non-inductively  wound  coil  of  manganin  alloy.   The 
pressure  response  of  this  alloy  is  extremely  linear,  and  can  be  cali- 
brated using  a  dead  weight  gauge  for  the  low  pressure  range  (where 
non-linearity  is  highest)  and  the  freezing  pressure  of  Mercury  at  O'^C, 
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which  is  7569  bars, to  fix  the  slope  in  the  linear  region.   The  gauge 
itself  consists  of  two  identical  coils,  one  at  system  pressure  and  one 
at   ambient   pressure,  set  in  a  constant  temperature  bath.   These  coils 
are  connected  to  a  Mueller  bridge  and  the  resistance  is  determined  by 
balancing  the  circuit  using  a  decade  box.   Again,  the  details  can  be 
found  in  Grindley  (1971). 

Unfortunately,  defects  in  the  decade  box  have  caused  the  manganin 
gauge  to  be  temporarily  abandoned,  and  replaced  with  a  Ruska  dead 
weight  gauge.   This  limits  the  pressure  range  to  12000  psia  (about  800 
bar) .   For  future  work,  the  decade  box  will  be  replaced  with  a  newer 
design  and  the  manganin  gauge  will  again  be  operational. 

Operating  Procedures 
As  in  the  previous  section,  the  operating  procedure  contains  little 
new  material  and  can  be  found  in  detail  in  Grindley  and  Lind  (1971)  and 
Grindley  (1971)  .   Only  the  change  to  the  use  of  the  dead  weight  gauge 
will  be  given  in  detail. 

Calibration  and  Correction  of  Raw  Data 

All  of  the  compression  and  tliermal  expansion  coefficients  for  the 

various  materials  are  readily  available  in  the  references,  but  the 

following  factors  cause  the  raw  data  to  need  correction: 

'  Thermal  expansion  of  lead  screw,  causing  the  measured  column 
length  to  be  too  short  at  higher  temperatures. 

'  Linear  compression  of  martensitic  plugs  and  floats,  causing  the 
liquid  to  have  more  volume  between  the  centers  of  these  pieces. 
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'Radial  compression  of  the  piezometer  tube  walls,  causing  the 
cross  sectional  area  of  the  fluid  column  to  increase  as  pressure 
increases . 

'Difference  in  pressure  at  the  pressure  measuring  point  and  the 
piezometer  tube  caused  by  the  height  of  the  mercury  transmitting 
fluid. 

Other  factors  influencing  the  measurements  are  automatically  accounted 
for  (e.g.  low  pressure  non-linearity  of  manganin  gauge,  and  temperature 
compensation  of  the  Heise  pressure  gauges — which  is  built  into  the 
gauges  themselves)  or  are  negligable  (e.g.  variation  in  bath  tempera- 
ture both  in  the  piezometer  bath  and  the  manganin  gauge  bath,  pressure 
effect  of  oil  column) . 

Actual  calibration  of  the  system  can  be  obtained  by  taking  measure- 
ments on  pure  water  and  comparing  it  to  available  data.   This  will  also 
lead  to  an  accurate  estimation  of  the  accuracy  of  the  volume  measure- 
ments. 

Piezometer  Preparation 

The  first  step  in  the  preparation  of  the  piezometer  apparatus  is 
the  placement  of  the  sample  into  the  piezometer  tube.   The  end  of  the 
tube  with  the  removeable  plug  is  pointed  upward,  the  plug  removed,  and 
the  liquid  slowly  injected  into  the  tube  using  a  syringe,  being  careful 
not  to  entrap  any  gas  in  the  tube.   The  float  is  then  inserted  and  the 
tube  filled  with  mercury,  again  using  a  syringe.   Reinsertion  of  the 
plug  follows,  and  mercury  should  flow  out  of  the  pressure  transmitting 
hole.   Additional  mercury  (about  1.5c)  is  placed  in  the  reservoir 
tube,  and  the  piezometer  tube  is  carefully  inverted  and  slid  into  this 
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tube.   The  entire  assembly  is  now  placed  into  the  pressure  tube,  taking 
care  to  push  the  reservoir  tube  as  far  into  the  pressure  tube  as  pos- 
sible.  The  pressure  tube  is  now  filled  with  oil  and  the  tubing  con- 
nections to  the  pressure  plant  reassembled. 

Pressure  Generation  and  Measurement 

The  sample  is  pressurized  by  operation  of  the  main  hand  pump  with 
the  valve  to  the  reversing  pump  closed.   Pumping  is  continued  until  the 
pressure  shown  on  the  Heise  gauge  is  approximately  the  desired  pressure 
(for  pressures  below  300  bar) .   The  pressure  can  be  adjusted  more  finely 
using  the  pressuring  adjusting  valve.   For  higher  pressures,  two 
changes  must  be  made.   Firstly,  the  pressure  must  be  generated  in  two 
steps,  and  secondly,  the  final  pressure  must  be  measured  with  the  raan- 
ganin  gauge  (the  procedure  for  the  dead-weight  gauge  is  given  later) . 

At  a  pressure  of  about  400  to  450  bar,  the  piston  of  the  pressure 
intensifier  reaches  its  limit  of  motion.   At  this  point,  the  valve  con- 
necting the  intensifier  to  the  rest  of  the  high  pressure  apparatus  is 
closed  and  the  valve  to  the  reversing  pump  is  opened.   The  valve  to 
the  oil  reservoir  for  the  main  hand  pump  is  opened  as  well.   The  re- 
versing pump  is  then  used  until  the  piston  has  returned  to  its  initial 
position  in  the  intensifier,  and  the  pressure  at  the  reversing  pump 
is  approximately  that  of  the  isolated  system.   The  reversing  pump  valve 
is  closed,  the  oil  reservoir  valve  opened,  and  then  the  valve  connecting 
the  intensifier  to  the  piezometer  apparatus  is  opened  very  slowly  to 
reduce  the  possibility  of  damage  caused  by  a  pressure  differential  across 
this  valve.   The  system  can  continue  to  be  pressurized  as  before. 
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Since  the  highest  range  on  the  Heise  gauges  ends  at  40,000  psig, 
the  manganin  gauge  must  now  be  relied  on  totally  for  pressure  measure- 
ments.  This  is  best  accomplished  by  setting  the  decade  box  to  the 
setting  which  corresponds  to  the  desired  pressure,  and  adjusting  the 
null  detector  to  a  low  sensitivity.   As  the  null  meter  shows  this 
pressure  being  approached,  the  sensitivity  is  increased  until  fine 
adjustment  is  to  be  made,  where  highest  sensitivity  must  be  used.   When 
the  meter  is  nulled  for  a  few  seconds,  the  pressure  tube  valve  should 
be  closed  to  reduce  the  chance  of  leakage,  and  the  volume  measurements 
can  be  taken. 

When  the  Ruska  dead-weight  gauge  is  used  instead  of  the  manganin 
gauge,  a  different  procedure  is  used.   In  this  case  the  desired  weight 
is  placed  on  the  gauge,  and  the  system  pressurized  until  the  weights 
are  lifted  well  above  the  reference  line  on  the  gauge.   The  weights  are 
allowed  to  settle  (due  to  slow  leaking  around  piston  in  dead  weight 
gauge)  and  the  pressure  tube  valve  is  closed  when  the  weights  cross 
the  reference  plane.   This  gives  the  desired  pressure  in  the  piezometer 
apparatus . 

Temperature  Measurement  and  Control 

The  temperature  of  the  system  is  controlled  by  tlie  bath  in  v>7hich 
the  piezometer  apparatus  is  situated.   This  bath  utilizes  a  silicon 
oil  which  has  a  very  low  vapor  pressure,  even  at  150"C.   The  temper- 
ature is  controlled  using  a  standard  proportional  controller,  and 
measured  using  a  Hewlett-Packard   digital  thermometer.   The  bath  is 
stirred  and  temperatures  vary  less  than  . 02' C  with  position  in  the 
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vicinity  of  the  pressure  tube.   Temperature  changes  can  be  accomplished 
quickly  using  two  hairpin  heaters  designed  for  that  purpose.   Because 
temperature  equilibrium  is  established  slowly,  the  entire  spectrum  of 
pressures  is  measured  at  a  given  temperature,  before  the  temperature 
is  changed.   Since  taking  measurements  for  the  entire  pressure  range 
is  a  day-long  procedure,  the  temperature  can  be  changed  and  allowed  to 
equilibrate  overnight.   In  any  case,  at  least  one  hour  should  be  allowed 
for  the  piezometer  apparatus  to  come  to  equilibrium  with  the  bath. 

Volume  Determination 

The  crux  of  this  experiment  is,  of  course,  determining  the  density, 
or  volume,  of  the  sample.   The  two  important  determinations  are  the 
length  of  the  liquid  column  and  the  size  of  the  sample.   The  composition 
of  the  sample  is  determined  when  the  sample  is  made  up;  however,  the 
weight,  or  number  of  moles  of  sample  used  is  not  directly  determined. 
This  is  because  some  solution  may  flow  past  the  float  during  the  filling 
procedure,  before  the  mercury  can  seal  the  interface.   The  mercury 
then  displaces  this  liquid,  which  spills  out  along  with  some  mercury 
when  the  plug  is  reinserted.   The  sample  size  then  must  be  determined 
by  comparison  of  low  pressure  measurements  taken  in  this  apparatus,  to 
any  available  low-pressure  data.   For  pressures  around  one  atmosphere, 
this  presents  little  difficulty. 

When  the  desired  temperature  and  pressure  have  been  attained 
inside  the  piezometer,  the  LVDT  is  used  to  measure   the  distance  be- 
tween the  upper  plug  and  the  float  of  the  piezometer  tube.   The  motor 
switch  is  turned  on  and  the  voltmeter  watched  for  the  presence  of  the 
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magnet  in  the  field  of  the  LVDT.   When  the  magnets  are  far  from  the 
coils,  there  is  zero  voltage.   In  the  vicinity  of  a  magnet,  the  voltage 
jumps  suddenly,  then  returns  to  zero  when  the  magnet  is  centered  in  the 
coil.   As  the  magnet  is  passed,  the  voltage  slowly  increases  and  then 
drops  sharply  after  the  magnet  passes  completely  out  of  the  coil.   The 
reading  on  the  turn  counter  for  the  motor  is  taken  for  the  upper  plug 
and  the  float,  and  is  converted  into  distance  using  the  gear  ratio 
from  the  motor  to  the  lathe  lead  screw,  and  the  screw  pitch  (corrected 
for  temperature).   The  tube  diameter  is  known,  and  after  it  is  con- 
nected for  pressure  effects,  can  be  multiplied  by  the  column  height  to 
give  the  volume  of  the  solution.   (The  measured  length  must  also  be 
connected  for  the  size  of  the  magnets,  since  it  locates  their  centers). 
Since  the  size  of  the  sample  can  be  determined,  the  molar  volume  and 
density  can  be  readily  calculated. 

Conclusions 
Although  the  experiment  has  been  operated,  no  real  data  has  been 
taken,  or  will  be  taken  by  the  author.   Only  the  construction  and  prep- 
aration  of  the  apparatus,  not  actual  measurements,  are  within  the 
scope  of  this  work. 


CHAPTER  7 
CONCLUSIONS 


This  work  has  focused  on  the  application  of  these  general  prin- 
ciples to  solution  thermodynamics  of  "reactive"  systems  — that  is,  the 
description  of  solutions  of  "reactive"  components  in  a  consistent  manner, 
using  both  theoretical  and  empirical  considerations.   First,  the  re- 
quirements of  a  consistent  model,  as  developed  from  classical  thermo- 
dynamics have  been  formulated,  both  in  the  form  of  a  test,  and  in  the 
form  of  a  procedure,  based  on  spatial  projections,  which  builds  this 
interval  consistency  into  the  model  development. 

Using  this  general  procedure,  the  solution  theory  of  Kirkwood 
and  Buff  (1951)  has  been  extended  to  these  "reactive"  solutions,  both 
in  terms  of  total  correlations  and  direct  correlations.   In  the  latter 
case,  a  set  of  first  order  partial  differential  equations  describing 
the  system  properties  in  terms  of  these  direct  correlation  functions 
and  reaction  extents  have  been  developed,  though  not  solved  explicitly. 
These  equations  rigorously  and  completely  describe  the  solution  pro- 
vided the  boundary  conditions  (reference  values  of  the  activities) 
and  equilibrium  constants  of  the  reactions  are  known. 

With  the  theoretical  basis  and  thermodynamic  consistency  firmly 
established,  these  have  been  applied  to  a  system  of  current  engineering 
interest — that  of  solutions  of  strong  electrolytes.   The  singular- 
ities caused  by  the  coulombic  potential  have  been  shown  to  be 
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rigorously  removed  by  the  projection  procedure.   Then,  a  composite 
model  based  on  the  primitive-ion  idea  has  been  developed,  but  only 
roughly  tested.   It  is  based  on  the  hard-sphere  equation  of  state  with 
the  limiting  laws  of  Debye  and  Hueckel.   The  form  originates  with  the 
mean  spherical  approximation  for  mixtures  of  ions.   Preliminary  cal- 
culations show  errors  of  less  than  five  per  cent  in  compressibilities, 
seven  per  cent  in  partial  molar  volumes,  and  eighteen  per  cent  in 
activity  coefficient  derivatives.   Improvements  are  expected  with  more 
careful  modeling  studies. 

Besides  these  calculations,  the  model  predicts  some  other  inter- 
esting results.   First,  the  Harned   coefficient  (to  the  lowest  order 
in  molality)  is  shown  to  be  a  weak  function  in  ionic  strength,  except 
under  some  restricted  conditions.   Secondly,  the  electrostatic  con- 
tributions from  Debye-Hueckel  and  the  mean  spherical  approximation, 
were  shown  to  be  expandable  in  a  power  series  in  Debye  lengths,  with 
the  parametric  form  of  the  first  several  terms  known.   The  dependence 
on  ion  charge  size  and  density  is  known  exactly;  only  the  coefficients 
are  model-dependent.   Thus,  coefficients  determined  from  data  on  one 
salt  can  be  used  for  all  salts.   These  can  be  compared  to  the  theo- 
retical values  obtained  here  for  the  mean  spherical  model.   Thirdly, 
this  formulation  takes  advantage  of  the  strength  of  the  fluctuation 
solution  theory.   The  derivatives  are  modeled  and  properties  are  ob- 
tained from  an  accessible  reference  state  and  the  change  due  to  com- 
position and  density  changes.   In  the  case  of  salts,  the  overall  concen- 
tration is  low  and  only  a  few  terms  in  the  concentration  expansion  need 
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to  be  retained.   It  is  expected  that  the  extension  of  this  formulation 
to  salting  effects  in  multicomponent  solvents  will  also  be  successful. 

Lastly,  an  experiment  has  been  built  and  a  procedure  established 
by  which  the  density  dependence  of  salt  solutions  can  be  determined. 
The  experiment  is  operational  and  procedure  specified.   The  method  for 
obtaining  activities  of  salts  at  high  pressure  from  these  and  available 
data  has  been  outlined.   This  information  allows  the  empirical  constants 
to  be  evaluated  from  the  experimental  data. 

In  all,  theoretical  basis  for  consistent  correlations  of  "reactive" 
solutions,  in  general,  and  for  those  based  on  direct  correlation  func- 
tions in  particular  has  been  established.   A  mathematical  model  which 
shows  promise  has  been  formulated,  and  an  experimental  apparatus  con- 
structed to  finalize  this  model.   But,  although  the  purpose  set  forth 
for  this  work  has  been  finished,  the  work  itself  is  more  of  a  starting 
point  for  the  engineering  applications.   There  are  several  directions 
to  be  taken  in  continuing  this  process.   First,  of  course,  is  the  deter- 
mination of  the  model  parameters,  and  testing,  or  even  revising,  it 
with  known  results.   When  complete,  an  explicit  model  for  compressibil- 
ities, partial  molar  volumes  and  activities  will  be  attained. 

From  here,  this  model  can  be  applied  to  various  systems,  an 
important  example  being  phase  equilibrium  of  salt-containing  solutions. 
Also,  this  model  is  applicable  to  systems  of  partially  dissociating 
salts,  but  in  this  case  the  differential  equations  of  Chapter  3 
describe  the  system  rather  than  the  more  simple  results  of  Chapter  4. 
A  solution  procedure  for  these  equations  must  be  developed. 
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A  totally  different  direction  can  be  taken  in  applying  the  totally- 
dissociating  component  equations  to  systems  other  than  strong  electro- 
lytes.  This  is  presently  being  done  for  the  "solution  of  groups"  model 
of  systems.   Of  course,  many  other  "reactive"  systems  could  be  described 
using  the  general  equations  of  Chapter  3,  provided  a  model  for  the 
direct  correlation  function  integrals  is  established,  and  a  general 
solution  procedure  devised.   Hopefully,  the  general  results  given  here 
will  indeed  prove  useful  in  predicting  properties  of  "reactive"  systems 
that  cannot  be  modelled  either  as  hard  spheres  or  primitive  ions. 


APPENDIX  A 
DERIVATIONS  OF  PRINCIPAL  RELATIONS 


Material  Balance — Equation  (2-1) 

The  component  vector,  N  ,  and  species  vector,  N^,  are  connected  by 

the  properties  of  a  set  of  reactions.   To  derive  the  expressions,  the 

first  group  contains  all  the  reactions  involving  only  the  components 

found  in  N   as  reactants.   This  labeled  the  first  set  of  reactions, 
— o 

with  r   reactions.   For  this  set,  N     is  defined  as  the  vector  of 
mole  numbers  of  all  n,  products  and  reactants  and  V    as  the  matrix 
of  stoichiometric  coefficients.   A  reference  component  is  chosen  for 
each  reaction  in  the  set  and  an  intensive  extent  of  reaction  matrix 
^    is  defined  by: 


r  ^(1) 


"  K-C 


extent  of  reaction  if  a  is  the  reference  com- 
ponent for  reaction  k.   This  quantity  indicates 
the  degree  to  which  the  reaction  has  proceeded. 
At  equilibrium,  its  value  is  found  from  the 
equilibrium  constant  of  the  reaction. 

if  a  is  not  the  reference  component  for  reac- 
tion k. 
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with  k  =  l,...,r  •  and  a  =  l,...,n  .   We  also  define  a  matrix  L    by: 


(L^^V  =  6.  =  s 
—   *ia    la 


1   i=a 


0   iM 


i-1, . . .n. 


a=l,  .  .  .n 


The  material  balance  for  the  first  set  of  reactions  is  then 


N^^^  =  (L^^)  +  V^^^C^^^N 


(A-1) 


.(1) 
Note  that  unlike  that  of  others,  the  present  ^    is  an  intensive 

quantity  and  its  form  depends  on  the  choice  of  reference  component 

for  each  reaction.   This  is  necessary  in  order  to  complete  the  desired 

manipulations  but  does  not  cause  any  difficulty  in  practice. 

Similarly,  we  have  for  the  second  set  of  reactions,  which  includes 

all  species  obtained  by  reactions  among  the  original  components  and  the 

species  created  in  reaction  set  1, 


(2)      (2)     (2)  (2)   (1)     (2)  ^   (2)  (2)    (1) 
N     =  (L    +  V   E,        )N    =  (L    +  V   g        )(L 


(1)  fl), 

+  V        E,        )N 
—   —    — o 


(A-2) 


where  the  matrices  superscripted  (2)  are  analogous  to  those  for  the 
first  set  except  they  are  of  dimension  n^xn-j^.   The  reference  component 
in  set  2  can  be  any  reactant  for  that  set,  including  a  product  of 
set  1. 


124 


Continuing  the  process  through  the  first  m  sets  which  constitute 
r  reactions  to  get  to  the  equilibrium  composition  N, 


,    (m)     (m)  rm)^^  (m-1)  _   (m-1)  (m-1)       (D 
N  =  (L    +  V   5    )(L      +  V     C     )...(L 


+  V   ^    )N   E  W  N  (A-3) 


Partial  Molar  Properties-Equations  (2-5)  and  (2-8) 
Starting  from  the  equations  of  equilibrium,  (2-2)  and  (2-4) 
we  write 


T     T       T     T 
dG  =  0  =  dN  y  +  N  dy  =  dN  y  +N  dy  (A-4) 

—  —    —   —      — O^Q   — O  -HD  ^  ' 


T        T 
where  N  dy  and  N  dy   are  zero  by  the  Gibbs-Duhem  equation.   Continuing, 
"  —     — o  ~o 

we  write  from  equation  (2-1) 


T       T 

0  =  A{\m  )^]^=  dNVy  +  N  d w  y  (a-5) 


Now 


T  (1)       (1)     (1)    T  (m-1)       (m-1)     (m-1)    T         (m)    T      (m)    T 

dW      =    (L        +V        C         )       • . . (L  +V  C  )     (dC         )     (V        )      + 


(1)       (1)     (1)    T  (£)    T      (£)    T  (m)       (m)     (m)    T 

+    (L        +V        5         )     ...(dC         )     (V        )     ...(L        +V        C         )      + 


(1)    T      (1)    T      (2)       (2)    (2)    T  (m)       (m)    (m)    T 

+    (dC         )     (V        )     (L        +V        ^         )       ...(L        +V        C         ) 
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So 


dw  y  =  0 


(A-6) 


The  above  equations  give 


T   T 

dN  (W  y) 

^o 


T 
dN  y 
— o^o 


Since  the  (dN  )  's  are  independent  at  any  N  ,T,P,  we  have  equation 

(2-5) .   Equation  (2-6)  is  obtained  by  noting  that  in  the  same  manner  in 

T 
which  dW  y  is  shown  to  be  zero,  it  can  be  shown  that 


(W^  y)  =  (J  y)  =  (y  ) 

—  —  a    —  —     — o  a 


(A-7) 


As  pointed  out  in  the  main  body,  the  n  -space  derivative  is  a  complex 
expression 


9N 


oa 


9N. 

=  y  — ^ 

i        oa 


3_ 
8N, 


T,P,N 


og^^a 


T,P,N 


o3?^a 


T,P,N.^^ 


=   1 


D      IP      Og 


9N 
i  oa 


9_, 
9Ni 


T,P,N 


oYT^a 


=  y  w. 

V     la 


9N, 


T,P,N. 


?^i 


T,P,N 


j^i 


9W 


+  yiN 


13 


^3     03   9N 


oa 


9_ 

9N, 


T,P,N       ,  T,P,N., 
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but  the  simple  expression  is  obtained 


and 


M   =  M        a  =  1,  . . . ,  n  (2-8) 

oa    a  '    '   o  ^   ^ 


Mixing  Properties-Equations  (2-12)  and  (2-18a) 
The  partial  molar  Gibbs  energies  of  mixing  can  be  evaluated  as 


AG^^^  =  y   -  y°  =   RTJlna  (A-8) 


AG          =   M      -  y°   =  RT!lna  (A-9) 

— o  ~o        — o  o 


so,   we  have 


T-— MIX           T             To                       T    o  o  o 
WAG         =Wy-Wy     =y     -Wy     +y     -y  (A-10) 
— o — o        — o 


which  is  equation  (2-12) .   The  general  form  of  this  equation  in  terms 
of  extents  of  reaction  and  equilibrium  constants  is  for  the  case  of 
equation  (2-14) 


T_^„X   —MIX      .   (1)  T    (1)      (1)   (1)  (1)  T   (2)  T    (2) 
W  Ag''^''  =  AG    +  R  T{(_^   )  £nK    +  (L   +v   £   )  (^   )  £nK    +. 

O 

(1)   (1)  (1)  T   (2)   (2)  (2)  T 
+  (L   +v   5    )  (L   +V   K        )       ... 


(m-1)   (m-1)  (m-l)  T   (m)  T    (m) 
x  (L     +v     K  )  (C    )  ^nK    }  (A-11) 


where  fcnK    is  the  vector  of  natural  logaritlims  of  the  equilibrium 
constants  for  the  reactions  in  set  £.   Although  appearing  complicated, 
this  is  a  straightforward  relation  to  use. 


APPENDIX  B 
EXAMPLES  OF  THE  PROJECTION  PROCEDURE 


A  General  Example 

We  first  illustrate  the  general  procedures  given  above  by  use 
of  an  example  which  should  indicate  the  types  of  manipulations 
commonly  encountered.   Then  we  describe  how  th  procedure  would  have 
been  implemented  in  a  case  from  the  literature. 

As  an  example  of  the  construction  of  W,  a  system  can  be  con- 
structed by  mixing  components  A  and  A  ,  which  undergo  the  following 
reactions: 


a)    A^  +  A^  ± 


2A, 


b)    A^  +  A^ 


=>   A, 


N   = 


ol 


o2 


Reaction  a)  makes  up  set  1,  and  reaction  b)  makes  up  set  2. 
The  matrices  for  set  1  are: 


,(1)  - 


-1 

-1 

2 
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(1) 

C   =  (C   0) 

—        a 


^     =   number  of  moles  of  1  reacted 
by  reaction  a)  divided  by  N 


ol 


(1) 


'l 

o' 

0 

1 

0 

0_ 

By  equation  (A-1) 


(1) 


'l 

0  " 

0 

1 

+ 

.0 

0  ^ 

25. 


^a 


2?    0  I 
*■  ^a     ■' 


1-e, 


<a        1 


'^l" 

.''02 , 

ol 


o2 


(B-1) 


For  the  second  set 


(2) 


(2) 


a^    0    0) 


C   =  number  of  moles  of  1  reacted 
by  reaction  b)  ,  divided  N^^-* 


(2) 


1 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

129 


(2) 


'1 

0 

0^ 

0 

1 

0 

0 

0 

1 

0 

0 

0 

-u 


0 

0  ' 

0 

0 

0 

0 

0 

0 

(!)■ 

N 

tl) 
^3 


? 

u 

1 

t 

0 
0 

0  1 
0 

1 

0 


(1)1 

3 


Combining  with  equation  (B-1)  the  composition  of  the  system  is  given 

in  terms  of   5  ,  C,  ,  N 
a    b  — o 


(2) 

N     =  N  = 


/■                 v 

N, 

1 

Nt 

2 

N^ 

= 

3 

N, 

I    4 

1-C- 


0 

r 


0 

1 

0 
0 


fl-^. 


2L 


0  • 

ol 

1 

0 

I    o2  J 

(1-?3)(1-V^^0l 

No2   -^aNol 
(2^a-?b+5a^b)Nol 

^b(l-5a)Nol 


(B-2) 


For  this  system  equation  (B-2)  gives  W  as 


W 


•(l-c,)(l-v 

0  ■ 

-^a 

1 

-^b(l-^'a>+2C, 

0 

.      ^(1-^a) 

0 

(B-3) 


Of  course,  W  would  be  different  if  other  reference  components  were 
used . 
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For  the  partial  molar  properties  equation  (2-8)  yields 


IL  = 


M. 


1  ~  9N 


9M 
8N 


ol 


T,P,N2,N3,N^ 


=  M 


ol 


T,P,N^2 


(B-4a) 


M^  =  M  „ 
^    o2 


(B-4b) 


From  equation  (2-9)  with  reaction  (a) 


M3  =  2  (^  +  M^)  =  i(M^^  +  M^^) 


(B-4c) 


From  equation  (2-9)  with  reaction  (b) 


"4  -  "1  +  "3  =  f"„i  *  ¥02 


(B-4d) 


For  mixing  properties  equation  (A-11)  with  equation  (2-15 


(2-15) 


gives 


(1-5,)(1-Sb) 


^a   -^b(l-?a)  +  2C3   C,(l-£  )  1 


b^-  -a- 
0 


f^l'''  1 
— ^MIX 

MIX 

AM. 

AM, 
^   ) 


ol 

'^a' 

1 

0 

0  ' 

K^^^^ 

+  R 

lO     . 

f(T£nK   ) 
a 

+   R  • 

I       o2 

.0 

1 

0  ^ 
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+ 


{-1  -1  -2} 


'E. 

0 

f(T£nK  ) 


(B-5a) 


-^^IX        — MIX 
leading  to  the  simple  result   A>C,    =  AM  ^    and 


o2 


AM^^^(l-5J(l-£  )  -CAm!,"^  -  {LiKJ    -   2C JAmV^"  +  Am;^"^^(1-5J 


'a'  "■      ^b   ^a  ^'2 


— MIX   — MIX 

,.v-  ■,  /    ^^  jAM,   +  AM,   4i,v-  ^.> 
b    a      a    3       4   b^   ^a^ 


— ^IX 


=  AM^^  +  R^^f(T£nK^)  +  RK^(l-E,^)f(T:lnK^) 


which  can  be  reduced  to 


—-MIX    T-MIX 
AM,     =  AM  , 
1        ol 


(B-5b) 


only  if  y    =  Pi •   The  simple  form  of  these  relations  occurs  because 
ol     L 

-—MIX      -r^IX 
of  cancellations.   However,  expressions  for   AM^    and   AM,     are 

complex  and  generally  not  of  interest. 

Finally,  for  excess  properties  and  activity  coefficients,  equation 

(2-20a)  becomes 


(i-^,)(l-Cb) 


b    a      aba 


gE' 

/"     _       ^ 

1 

c^ 

— 

ol 

2 

G^ 

gE 

3 

I    o2 

C^' 

^^■4   J 
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+   R   T 


2.nK      +  RT  \ 
a 


1 

0 

0" 

+ 

<^.' 

0 

1 

0 

[0    J 

{-1      -1      2} 


''=1 

0 

10     J 

InK, 


f  i\  \ 


£.nx 


+  R  T 


ol 


£nx 


o2 


fd-yCl-C,)      <^      -^b^l-V    +   ^C,      ^b^^-^a^ 


Jinx 


?,nx 


£nx 


After  considerable  algebra,  this  reduces  to  what  would  be  expected 


G^  =  G^  +  R  T{£nx   -  £nx^}    a  =  1,2 
ot     oa    g     oa      c 


(B-6) 


For  the  reference  states  given  above 


£nY   =  £nY  +  £n(x  /x   )     a  =1,2 
oa     a      a  oa 


with 


X  = 
1 


M 


I  W,  N 

la  oa 
a 


yy  w.  N 

V^  la 
la 

1  -  x^ 


X^T  E  N  i/y  N 
ol    ol  ^   oa 


oa 


a 
=  1  -  X 


o2 


Then,  for  this  example 


ilnYg;!^  ^  ^'"^1  ~  ^^ 


l-^b(^-^a)-ol 

Id-c^xi-^b)   J 


(B-7) 
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If  the  extents  of  reaction  are  arbitrarily  set,  equation  (B-7)  can 
be  used  directly  in  equation  (2-21a)  along  with  a  model  for  J^ny,  •   If 
"equilibrium"  is  assumed,  then  Cg'  ^u'  ^l»  ^2'  ^3  ^^^   ■^A  must  be  found 
from  equations  (2-1)  and  (A-11)  specified  values  of  x  ■.  ,  K  and  K,  ,  and 
a  model  for  the  Y-;  •   In  this  case,  the  expressions  are 

^     °^  (B-8a) 


^'  ^-h^'<aKl 


l-(l+^a)>^o1 

^2  =  1-(1-?  )C  X  ^  (2-Sb) 

a  b  ol 

{2g^  -  Cb(l-ga)Kl 
■^       ^  ^a   %   ol 


X,  =  1  -  X-]  -  X2  -  x„  (B-8d) 

2 
Kg  =  (Y3X3)  /(Y^x^Y2X2)  (B-8e) 


^b  =  V4/^^1^^3"3^  ^^"^^^ 


An  Example  of  Continuous  Association  in  the  Chemical  Theory 

Prausnitz  (1969:  Sections  7.13  and  7.14)  gives  two  examples  of 

the  chemical  theory.   In  the  first,  all  degrees  of  association  are 

allowed  among  an  alcohol  (2)  in  a  hydrocarbon  (1).   Thus,  starting 

with  N  ,  moles  of  hydrocarbon  and  N  „  moles  of  alcohol,  there  are 
ol  •'  o2 

continued  successive  reactions 


^  +  \  <=='  Vl   "  =  2,3,...,^ 
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If,  for  each  reaction  (£) ,  we  choose  the  reference  component  to  be 

(I) 

the  n-mer,  then  the  mole  vectors  have  elements  (N   )  for  n>2,  t>l. 

—   n      —    — 

The  stoichiometric  vectors  have  elements 


(V<»,,  =  -2   ;  (V<»,3  =  1 


(I). 


u).         ,  a). 


(V   )2  =  -1  =  (V   )^^i  ;   (V   )_^^2  =  ^  ' 


(l) 
(V   ).  =  0  ,  i?^2,  £+1,  1+2;   l>2 


The  extent  of  reaction  vectors  have  elements 


(n-1) 

a         ) 


=  K 


,(n-l) 

^n+1 


n-1  -  „(n-2) 


(i^^^k  =  °'   ''^  ^"^^'  "^^     ^'^"^^ 


n 


and  two  important  stoichiometric  relationships  are 


N. 


.(£) 


1  -C. 


J?,    k 


2+    I    T  r' 

k=2  k=2 


■  N  T    S,>2       (B-lOa) 
o2    — 


(£) 
V+1 


lk=l 


o2 


£>1       (B-lOb) 


The  projector  matrix  is  then 


T 
W   = 


1     0 


N. 


o2 


(1  -  ^2)^1 


(l-^^-l>JJ   ^k 


k=l 


(B-11) 
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where  N.  is  the  value   lim  (N   )  .  =  {l  -  ^,  (2  +  E   IT  ^,  , }   N  ^ 

£-  -    "         ^    k=2  k'=2^     °2 


While  these  could  be  used  in  equation  (2-16a) ,  here  their  use  is  much 
simpler  in  equation  (2-16b)  and  thus  in  equation  (2-22). 

The  expression  chosen  for  the  Gibbs  free  energy  of  mixing  of  the 
species  is  obtained  from  two  contributions  (Renon  and  Prausnitz,  1967) 


^^IX^   ^-MIX    ^   ^-MIX 


The  first  (chemical)  is  for  forming  a  lattice  of  hydrogen-bonded 
alcohol  species  at  the  volume  fraction  of  the  solution  from  a  reference 
state  of  the  pure,  crystalline,  oriented  species.   The  second  term 
(physical)  is  a  simple  Van  Laar  expression  for  interactions  on  the 
lattice  in  the  solution  of  species.   These  yield  for  the  species 
(including  hydrocarbon) 


(AG^^^)/RT  =   {(AG^^^.  +  (AG™)^}/RT 
—    1  — c    J-    — P    i 


{£n$  +  (1-0  )  -  V   );   N./  I     V.N,} 

\=2  j  =  l   -J  J 


1=2       1=2 


oo     oo 


y    y  B.A.^.}/RT  (B-12a) 

i=2  j=2  ^J  ^  J 
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(AG^^^)  /RT=   {(AG^^^).  +   (AG^I^).}/RT 
—    i        —c   1       P   1 

V.  y  N. 

v^                      2=2 
=   In  ($./v.)   -  — {$,  +  ln(T+l)  -1}  -  — 

^   ^      '1   '  I  v., 


.N. 
J  J 


j=l 


+   ln(T+l)  +  In  (2v  ) 


+  V.  {  $,(6^  .  -   Z   3^  .$.)  +  Z   $.(26.. 
1    1   li   ^^2   IJ  J     j=2   J    ^J 


Z  3  ,$.,)} 

j'=2   JJ   J 


(B-12b) 


Using  equation  (2-22) ,  the  desired  activity  coefficients  are 


V  ,  <» 

In  Y  ,  =  ln($  /x   )  +  (!-$)-  Z  <D  /(v./v  J 

ol       1  ol  1    ^  0  .  T  1   1   o2 

oz  1=2 


'ol 


CO       CO 


+  -^   {(1  -  $  )   E   g  $.  -  I   E   B.  $.<!>.} 
^^  i   .  _   ll  1    .  „  .  „   li  1  1 

1=2  j=2   -J    •> 


i=2 


(B~13a) 


In  Y  ,  =  In  (<I'„/<f°x   )  -  -^  $   -   Z   ($.  -  $°)/(v./v   ) 
o2        2   2  o2    v^^   1    j=2   J     J     J   °2 


j=2 


RT   ^^l^Pl2    .-.  -Ij  y 
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+  E  [  2B^.($  -  $°)  -  E     («>.$. ,  -  $°$°')  ^,,^^}  (B-13b) 


Other  than  the  basic  models,  no  assumptions  have  been  made.   Yet  the 

results  are  readily  found  from  (2-22)  . 

Major  simplifications  appear  if  it  is  assumed  that  the  volumes 

are  additive 

v./v  ^  =  i  -  1  (B-14) 

1   o2 

that  the  physical  interactions  are  independent  of  polymerization 

^-  -  1^12'   J-^  '   ^ij  -  ^22  '  ^'J-2  ^^~^^^ 

and  that  the  polymerization  reactions  are  independent  of  polymer- 
ization and  are  expressed  by 


Cj^  =  K$2  =  ^    ;  k  >  2  (B-16) 


where  K  is  a  function  of  temperature  only.   (Because  of  the  difference 
of  the  dimerization  reaction  E,^    is  different  from  ^.)   These  assumptions 
used  in  equation  (B-lOb)  yield 


$  =  (i-T)(K$^)^  $^        i  >  3  (B-17a) 

i  2      2 


CO  fy 

$   =   E   $.  =  <J)^/(1  -  K<J)„)  (B-17b) 

o2   i=2   ^    2        2 
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E   $./(i-l)  =  ^J    (1-K<D  )  (B-17c) 

1=2   ^  ^       ^ 


A  value  for  $„  can  be  found  in  terms  of  K  and  0  ^  from  (B-17b) .  A 
■^  o2 

value  for  ^      can  be  found  from  equations  (B-17b)  and  (B-lOa) ;  here  it 
is  K<I>o(l-K$«).   Using  equations  (B-17)  in  (B-13)  yields  the  results 
of  Renon  and  Prausnitz  (1967) .   The  difference  in  the  derivations  are 
that  we  do  not  need  to  define  the  relationship  for  the  equilibrium 
constant  K  in  terms  of  any  model  and  that  the  procedure  is  independent 
of  the  particular  lattice  model  chosen.   Introducing  alternative  assump- 
tions about  the  reactions  (e.g.,  1st  K  for  the  dimerization  reaction  be 
different  from  that  for  other  reactions)  would  change  only  equations 
(B-17)  and  the  precise  expressions  are  easily  found  from  equations 
(B-13). 


APPENDIX  C 
EXPRESSIONS  FOR  THE  UNIQUAC  AND  UNIFAC  MODELS 


In  both  the  UNIQUAC  and  UNIFAC  models,  the  excess  Gibbs  energy  is 
the  sum  of  two  terms. 


o     o  combinatorial     o  residual 


(C-1) 


For  both  models,  the  combinatorial  parts  are  the  same.   For  UNIQUAC,  the 
parameters  r   and  o   are  for  molecules  while  for  UNIFAC,  they  are  ob- 
tained by  projection  of  the  group  parameters  R.  and  Q.. 


q  =  J  V .  Q. 
VL  ^   ia^i 


(C-2) 


(G^/RT) 


.  -,  =  -  (S^/R)   ,     ■,  =  y  N  £n((j)  /x) 
combinatorial  athermal   ^   oa   ^a   oa 

a 


In  UNIQUAC 


(^o/^^)  residual  =  "  I    ^ca^^'^^j   %    exp(-Au3^/RT)  ] 


(C-3b) 


a 


and  <t)      and  6   are  volume  and  surface  fractions  for  the  component  a, 

a     a 

In  UNIFAC 


o 


,RT, 


mk 


residual     a    ^k      "■   m 
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(C-4) 


(a) 
where  0  is  the  group  surface  fraction,  0    is  0  in  a  solution  of 
m  mm 

pure  a. 

In  the  UNIQUAC  model  for  molecules,  the  total  interaction  energy, 
U,  is  given  by 


u  =  ^zyN  q   ye„   u„ 

2   ^  oa  a  ^  3a  Ba 
a      p 


where 


6a 


N^gqgexp(-Aug^/RT)     N^gg^T^^ 


(C-5) 


Au   =  —  [u   -  u   ] 
ya   2   ya    aa 


In  the  UNIFAC  model  for  groups,  U  is 


U  =  4  Z    y    N.Q.    \   0. .U.. 
2        h     11    .      J 1  ji 


(C-6) 


with 


N.Q.    exp(-AU. ./RT)  N.Q.T. 


0. 


J_J 


JA 


Ji 


■1    .1    Ji 

I   N^Q^   exp(-AU^./RT)         I   N^Q^T^. 
k  k 


AU,  .    =  -   [U,  .    -   U    .] 
ki        2        ki  11 


The  partial  molar  energies  for  UNIQUAC  are 


q  N  „q„u„  exp(-Auo  /RT)  ,      q  N   q  u   exp(-Au^^/RT) 

-  ^  \      y  ^n  oB  B  Bn     Bn  j^\_  y  ^n  og^g  na     na 

°'^    ^    B    y  N   q  exp(-Au   /RT)  ^  a    I  N   q  exp(-Au  ^^/RT) 
'^     oy  Y      YH                 oY  Y      Ytt 
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-—   zYqN      q      y 
a  3 


'I   N^^q^  exp(-Au^^/RT)l 


(C-7) 


and    for   UNIFAC 


^        ^        j         ^  N  Q   exp(-AU      /RT)         ^        i        ^  N  Q      exp(-AU.  . /RT) 


'k^r 


kr 


k^k 


ki' 


N.Q.U..    exp(-AU. ./RT)    exp(-AU. . /RT) 

+  7^1  Q.N.Q.  I  ^LJ_Ji_ _J^ Al 

•  fl   QA   exp(-AU^./RT)]' 

k 


(C-8) 


To  test  the  projections  of  equation  (2-32),  we  must  first  make  the  stan- 
dard proiections,  q   =  T  V   Q   and  N„  =  T  v„  N    to  write: 

p  ^  '^        e 


\  <' 


^o^^n^-    N  ^Q.U.„  exp(-AU.„/RT) 


=  i.  7  yyy  i?-n  ^  .jg  og^j  jg, 


ia: 


^aj   II  \/oy\   exp(-AUj^^/RT) 


V„  Q„V.  N  Q.U„.  exp(-AU../RT) 


^j    II   ^wN.v^V  exp(-AU^./RT) 
ky 


'kY'"oY~^k  """''  ^  "  k j  ' 


(C-9) 


V„  QnV.  N  Q.v.qN  QQ.U..exp(-AU../RT)exp(-AU   /RT) 


£gjBi 


'yy  V,  N  Q,  exp(-AU,  ./RT)" 
^^  ky  oy  k   ^    ki 


while 


on   2 


v„  Q„v.  N  Q.u   exp(-Au  /RT) 

yyy  ill  i?^  .ig  oa  J  an  ^         an 

£gi    yy  V,  N   Q,  exp(-Au   /RT) 
•J    i^   I<.y  oy^k   '^    yn 


v„  Q„V.  N  Q.u   exp(-Au   /RT) 
1   yyy  g.n^  .la  oa  J  na       ng 


(C-10) 


^ctj    II   ^,  N  Q,  exp(-Au  ^/RT) 

■J        LL     i^y  oy  k       ya 
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\  -  mil 

£aj3i 


V.    Q„V.    N     Q.V.qN   oQ.u„      exp(-Auo    /RT)    exp(-Au      /RT) 


'yy   V,     N      Q,     exp(-Au      /RT) 
f:^     ky  oy  k       ^         ya 


It  it  apparent  that  the  two  equations  cannot  be  made  identical, 

because  the  AU...  etc.,  will  not  be  the  same  as  Au   ,  etc.,  in  the 
xj      '  an 

summations.  If  the  projections  of  q  and  N  are  made,  the  "Boltzmann- 
factor  weighted"  projection  of  the  surface  fractions  cannot  be  made, 
without  a  loss  of  consistency. 


APPENDIX  D 

CONSTRUCTION  OF  MATRICES  FOR 

GENERAL  MATERIAL  BALANCE 


Here  we  consider  the  material  balance  and  associated  matrices 
in  the  general  case  of  subsequent  reactions.   First 


N  =  WN 
-   — o 


where 


[l(2)h.v(2)^(2)jjJ1)^  ^(1)^(1)^  (,_,) 

as  derived  in  Appendix  A.   When  the  matrices  are  divided  into  blocks 

as  defined  in  equations  (3-3)  to  (3-5),  equation   (D-1)  can  be  rewritten 


where 


W 


U) 


I  +  v^^'h^'^ 


U)  (a) 

-D  - 


^0)^H) 


n  X  n„  - 


(D-2) 


(£)         (£) 

By  analogy  to  equations  (3-8)  and  (3-13)  we  define  Z    and  Y 
by  the  equations 


1A3 
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,(0 
-o 


T        T 

,  ii)\-l     (JO 


"£  ^  ^"r^ii-r'^^     ^°-3^ 


rW    _ 


T        T 
,  (£)\-l   (ii)^ 

0 


"S,  ^  ^£-1 


(D-4) 


As  before, 


and 


-o   -      - 


y(0\(£)  ^  J 


(D-5a) 


(D-5b) 


Now  the  overall  Y  and  Z  matrices  can  be  constructed.   The  Y  matrix  is 
-     -o 

the  more  straightforward,  being 


-  o 


(D-6) 


From  equation   (D-5b)  it  can  be  easily  seen  that  Y  W  =  I  as  required, 

The  construction  of  the  Z  matrix  is  less  obvious.   The  first 

-o 

/         \    T         r,(m) 
n  X  (n-n   ^-r  )  columns  are  Z    since 
m-1  m  -o 

The  next  n  x  (n  ^ -n   _-r  ,)  columns  of  Z  are  y^^^z^"""'"^ 
m-1   ra-2  ra-1  -o     -   -o 


^^(m-1)  ^(m-iym-2)___^(l)  ^^ 
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where  equations  (D-5)  have  been  employed.  Continuing  in  an  analogous 
manner,  the  Z  matrix  can  be  constructed  by  juxtaposing  m  submatrices 
and  has  the  form 


-o    -o   '  -   -o     »•••»_   _     •••_     _o   '  ' 

^im)^(.m-l)  ^    ^^(3)^  (2),  Y^Ydn-l)    y^^^^^^Z  ^^h 


(D-7) 


The  total  number  of  columns  of  Z  is 

-o 


(n-n  1-r  )  +  (n   -n   -r   )  +  ...  +  (n  -n   -r  )  + 

m-1  m      m-1  m-2  m-i  I      ic-l 

m 


(n^-n^-r^)  +  (n^-%-r^)    =  n-n^+  J 


r  =  n-n  -r 

£=1  ™ 

So,  the  dimensions  of  Z  are  n  x  (n-n  -r)  as  required. 

-o  o 

The  U  and  Z^  matrices  cannot  be  written  explicitly  in  the  general 
case  and  so  are  simply  defined  by  equations  (3-9)  and  (3-11).   However, 
they  are  not  explicitly  used  either. 


APPENDIX  E 
DERIVATIVES  WITH  RESPECT  TO  THE  COMPONENT 
CHEMICAL  POTENTIALS 


The  derivatives  of  the  grand  partition  function  (5)  and  the 
average  composition  vectors  (<N>  or  N  )  with  respect  to  component  chem- 
ical potentials  (u  or  InX)    are  treated  in  detail  below.   The  difficulty 
in  evaluating  these  derivatives  stems  from  the  fact  that  the  vector  of 
Lagrange  multipliers  in  the  constraint  space,  g.n9,  is  by  necessity  a 
function  of  those  in  component  space,  InX. 

Allowing  V  to  represent  the  matrix  of  stoichiometric  coefficients 
of  all  n  species  in  all  r  reactions  (regardless  of  reaction  set) ,  the 
equations  of  reaction  equilibria  are 


T 
V  M  =  0 


(2-4) 


or,  partitioning  the  matrices  as  in  equation  (3-5) 


,  T  T   T. 
(^C  ^D  Hl^ 


He 


Hd 


T    ,   T    ,   T 
=  ^cHc  +  Vd  +  ^iHl  =  2 


(E-1) 


Differentiating  with  respect  to  £nA   (=3y   )  and  noting  that  u^  =  u  as 

a     oa  °      -C   -D 

per  equation  (2-6) ,  equation  (E-1)  becomes 


T     ,   T 
-C  -a        -D 


■^Hd 


dlnX 


+  ^I 


^Hi 


dinX 


T,V,X 


=  0 


(E-2) 


y^a 


T,V,X 


Y?^a 
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Using  the  fact  that  v   is  nonsingular  by  construction,  we  premultiply 


T  -1 
by  (v_)   and  rearrange  to 


^Hd 


9  InX 


T,V,X 


,  T--1  T     ,  T,-l  T  ^^-I 

c 


(E-3) 


T,V,A 


y^a 


Next,  equation  (3-25)  is  differentiated  with  respect  to  3  ^nX  , 
and  again,  the  y  vector  is  partitioned,  giving 


8m 


dlnX 


T,V,X 


y^oL 


e 

-a 


^Hd 


8JlnA 


a 


^Hi 


T,V,A 


y^a 


dln\ 


T,V,A 


Y?^a 


(E-4a) 


8U 


8£nA 


9Z 


toA  +  (U)   + 


-  a    92,nA 


T,V,A 


a 


9£ne 


?,nA  +  Z 


-  92,nA 


y^a 


T,V,A 


a 


Y^a 


T,V,A 


(E-4b) 


Y^a 


Substituting  equation  (E-3)  into  (E-4a)  gives,  after  rearrangement. 


8y 


9S,nA 


T,V,A 


Yi^ct 


I 

.   T.-l   T 

e     + 

-a 

0 

^ 

.  T--1  T 


^Hl 


3  HnX 


a 


T,V,A 

Y 


5^a 


9^1- 


=  (Y)   +  Z   3 


-a   -o    9£nA 


(E-5) 


T,V,A 


y^oL 
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Finally,  combination  of  equations  (E-4b)  and  (E-5)  leads  to  the  depen- 


dence of  J?,n6  on  the  JlnA  's,  namely 


9g,n9| 

3£nX 
a' 


=  (Y)„,  +  Z 


^Hi 


-  ct   -o   9£nX 


T,V,A 


y^a 


a 


T,V,A 


9Z 


dlnX 


9U 


Zne   - 


dlnX 


HnX 


T,V,X 


y^a 


T,V,A 


y^oL 


(E-6) 


where  the  notation  (Y)^  or  (U)^  indicates  the  vector  which  comprises 
the  ct-th  column  of  Y  or  U. 

To  examine  the  fluctuations,  we  look  at  the  derivatives  of  the 
grand  partition  function.   First,  differentiation  is  carried  out  with 
respect  to  £nA 


1      9H 


5    9JinA 


■  i  H  h' 


a' 


T,V,A 


y^a 


N   q 


9U 


dUnX 


9Z 
lEl  +   (U)     +  ^-— 


ind 


T,V,A 


y^a 


T,V,A 


y^o. 


9£n8 


+   Z 


-    9£nA 


T,V,A 


Y?^c 


exp[-eE^^   +  N    (U    taA      +  Z    £n9)]  (£-7) 


which, with  equation  (E-6),  reduces  to 


1      9E 


H    9JlnA 


=   <N> 


T,V,A 


Y^a 


^Hi 


(Y)      +  Z     ^ , 
-  a       -o      9jlnA 


T,V,A 


Y^a 


=    <N>  (Y)      =   N 

-        -  a  oa 


(E-8) 


U9 


where  equations  (3-7)  and  (3-14a)  were  used  to  obtain  the  last  two 

equalities.   This  result  is  what  is  expected  since  inX      =  gu 

a     oa 

The  quantities  of  greatest  interest  are  the  derivatives  of  average 
composition  with  respect  to  component  chemical  potentials. 


8<'N> 


9JlnA 


a^nX 


T,V,A 


Y?^a 


T,V,X 


{  4  5^  );  N  exp[-eE   +N^(y  £nA  +  Z  RnQ)]  } 
"  N  q  -^ 


y^a 


<nn'^> 


3U 


3£nX 


9Z 


inX   +  (U)   +  ^ 

-  a   9 


2.nX 


£n9 


T,V,A 


Y/a 


T,V,X 


y^a 


9£ne 


+  Z 


-  9£nX 


T,V,X 


Y?^a 


-  <N> 


1  9; 


dlnX 


T,V,X 


y^ct 


(E-9) 


Equations  (3-26) ,  (D-6)  and  (D-8)  allow  this  to  be  simplified  to 


9<N> 


9K,nX 


9<N> 


96p 


T,V,X 


oa 


T 

<NN  5 


Y?^a 


T.V.u, 


Y?^ct 


^Hi 


(Y)   +  Z  - 

-  a   -o  9  p 


oa 


T,V,y^ 


3   , 

Yf« 


-  <N^ 


oa 


(E-10) 


Finally,  equation  (E-10)  can  be  projected  into  component  space  using 
equation  (3-14a)  and  the  fact  that  Y  is  a  constant  matrix,  giving 


9<N> 


9Bm 


oa 


9y'^<N> 


9Bii 


T,V,p^ 


oa 


9N 
-o 


9  3  U 


Y^cx 


T,V,M, 


oa 


C^a 


T,V,M^ 


jT^a 


150 


which  becomes  equation  (3-36),  and,  since  <N>  Z  =0,  yields  equation 

(3-37). 

T 
Another  important  relation  that  arises  from  <N>  Z  =0  and  the  fact 

-  -o   - 

that   Z      is   constant   is 
-o 


8<N> 


-o   3  By 


oa 


T,V,y, 


=  0  =  Z     — ^; — 
-o      <N> 


y^oL 


^Hi 


(Y)    +  Z      ^ 

-tt       -o    9y 


oa 


T,V,y 


Y?^ct 


(E-11) 


This  relation  is  used  in  Chapter  4, 


APPENDIX  F 

PROJECTION  OF  INVERSE  MOLE  FRACTIONS  OF  LIMITING 

REACTANTS  IN  THE  LIMIT  OF  COMPLETE  REACTIONS 


Since  x7  diverges  in  the  limit  of  complete  reactions,  it  is 
necessary  to  first  project  x7  by  (W+K)  at  non-zero  limiting  reactant 
concentrations,  then  perform  the  limiting  process.   This  is,  of  course, 
the  rigorous  way  to  evaluate  the  limit,  but  was  unnecessary  for  the 
remainder  of  (X  -C) . 

In  the  general  case  (c's  unspecified), 


(W+K)^ 


n  n     n 

o  o     o       9  5 

^       a=n  +1        a   y=1   a=n  +i      of 


y^. 


T,V,N, 


yH 


n 


=  '^i6''  ^   ./ia 


a=n  +1 
s 


3  5. 


a      a  of 


T,V,N 


V3 


(F-1) 


and 


0 


9  h' 


0    0 
0 
0    0    0 


IJ 


i  6.. 


if  n  +1  :2  i  <  n 
s         o 


otherwise 


N6 


IJ 


N  .+v. .f .N 
oi   11^1  or, 


n  +1  <  i  <  n 
s         o 


otherwise 


(F-2) 
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The  projection  of  the  divergent  matrix  yields 


(W+K) 


0 

0 

0 

0 

h' 

0 

0 

0 

0 

(W+K) 


=  N   I 


i=n  +1 
s 


n 


5-0+  I  V.  (?  5„ 

i6      , ,   la   a  Br 
a=n  +1  a 

s 


By 


H, 


+  N 


or   8N  , 
a    o£ 


T,V,N 


OY^B 


ir    ,  ^    la'   a   Br  , 
a  =n  +1  a 

s 


+  N 


^^a 

"'^a- 

'\b 

T,V,N 


(N  .+v. .C.N   ) 
oi  11  1  or. 


OYT^B 


(F-3) 


As  ?   approaches  its  value  for  a  complete  reaction  (equation  3-63) 


-       1    oa 

5  = t; e 

a     V  N       a 
aa   or 


(e   ^  0) 


(F-4) 


where  e   is  positive  and  goes  to  zero  as  complete  reaction  is  attained. 


Also 
9C 


8N 


oB 


N   6  „-N   6„ 
,    or  aB   oa  Br    3^ 
la         a     a 


T,V,N 


aa      N 


3N 


oYT^B 


or 


oB 


(F-5) 


T,V,N 


OY5^6 


1   N  3C 

I    oct  ct 

Since  E.      is  bounded  (0  1  C   < t; )  and  the  limit  of  t-;— 

a                 a     V    N  3N 


aa   or 


oB 


T,V,N 


OY?^B 


3?, 


exists,  the  quantity 


3N 


oB 


must  also  approach  zero  as 


T,V,N 


OYT^B 
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e  approaches  zero.   So,  for  small  enough  e 


8e 


8N 


oB 


=  M„  e  /N 
pa  a  or 


(F-6) 


T,V,N, 


a 


Y^B 


3e 

where  M„   is  always  finite  and  of  the  same  sign  as  rrr 

Ba  3N 


oB 


T,V,N, 


YT^a 


Substitution  of  the  results  of  equations  (F-4),  (F-5)  and  (F-6) 
into  equation  (F-3)  yields  (for  E,    's  very  near  their  complete  reaction 
values) 


(W+K) 


fo   0   o1 

0    x;;^    0  (w+'^) 

0    0    0 

J  J 


n 
o 

i=n  +1 
s 


(f 


iB   ^   ,-,    V   aB  la  a  gr 
a=n  +1    aa  c 

s 


3y 


)  f 


-  V.  M  e  ) 

la  Boi  a 


6.  +  y    ( 5  ,  -V.  ,e  ,6  ,  -V.  ,M   ,c  ,) 

XY      ,   V  ,  .ay   la  a   a  Y   la  yc   c 
a=n  +1   a  a 


ii» 


(N  . ^  .-e.N    )  =  N  y 

oi  V. .  oi  X  or .      .  ^   , , 
XX         1     x=n  +1 

s 


6.  -  -^.  -v..£.(6-  -M^.) 
x3  V. .  xB  11  1  3r .   Bi 


11 


X  [6.  --^6.  -v..e.(6   -M  .)]  v(-e.N   ) 


XY  V .  .  XY   XX  X  yr .   yx 

XI  1 


X  or . 

1 


"o    \i(^Br.-V(\r.-V 

=  -Ny     i-.^ i c, 


i=n  +1 
s 


N 


(F-7) 


or . 
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To  take  the  limit  of  complete  reactions,  the  limit  as  all  the  e  's  go 
to  zero  is  taken  on  the  right  hand  side  of  equation  (F-7) .   This  leads 
to  the  result  that 


(W+K) 


0 

0 

0 

0 

h' 

0 

0 

0 

0 

(W+K)  =  0 


(F-8) 


in  the  limit  of  complete  reactions. 


APPENDIX  G 

ALTERNATIVE  FORMULATION  OF  FLUCTUATION  PROPERTIES  FROM 

INTEGRALS  OF  THE  RADIAL  DISTRIBUTION  FUNCTION 


Constraints  on  the  Radial  Distribution  Function  Integrals 
Friedman  and  Ramanathan  (1970)  place  the  following  constraints  on 
the  integrals  of  the  radial  distribution  function,  based  on  charge  neu- 
trality and  material  balance  arguments: 

(z).,  +y  (H)..x.(z).,  =0     l^i%,   1-k-n-n      (G-la) 
-  xk   ^  -  13  J  -  jk  o 

or 

(  2  +  H  X  )  _z  =  0  (G-lb) 

There  are  two  consequences  of  the  presence  of  these  constraints.   The 
form  of  the  constraints  given  in  equation  (G-lb)  makes  it  clear  that  the 
relationship  between  A_   and  ]I  differs  considerably  from  equation  (4-23) , 
in  fact,  it  becomes  a  much  simpler  relationship.   The  form  of  equation 
(G-la)  shows  that  the  models  used  for  the  (H) . .  must  be  related  to  each 
other  in  this  fashion.   This  implies  that  the  constraints  on  the  (H) . . 
may  have  to  vary  from  system  to  system, particulary  when  common  ions  are 
involved.   Examples  of  both  of  these  effects  are  given  below. 

Equation  (G-lb)  can  be  used  with  equation  (4-21)  to  give  the  fol- 
lowing relationship  between  A_   and  H: 

A_  =  y^    (X  +  XHX)  X  (^"2) 

In  terms  of  the  desired  quantities,  this  becomes 
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A      =[Y^(X  +  XUX)y_]^ 


(G-3) 


Although  this  is  much  simpler  than  equation  (4-23)  and  its  inverse,  the 
need  of  a  matrix  inversion  in  equation  (G-3)  still  makes  determination 
of  the  elements  of  A   much  more  difficult  in  terms  of  the  (H) . .  than  the 
(C^  ) . . ,  and  this  difficulty  increases  greatly  as  the  number  of  compo- 
nents increase. 

Single-Salt,  Single-Solvent  System 
In  this  simplest  case,  A_   becomes 


A   = 


10    0 
0   l/v_^  0 


^ol+^ol  ®  11   ^ol^+®  1+   ^ol'^-  ^^)  1- 


^ol^+®l+ 


X  ^x  (H)^ 
01  -  —  i- 


V  ^+®- 


^+^_(H)4- 


x^x_(H) 


X  +  X  (H) 


0    0 
0    1/v^ 
0    0 


^01+  ^01®  11 


\l\2^^h+ 


\l\2^^h+ 


o2    2  ,  , 


(G-4) 


Equation  (G-la)  shows  that 


(H)^_   = 

^^\+ 

fu^ 

1 

<■!)+_  - 

V     X     „ 

-   o2 

+  (H) 


V  X 

+  o2 


+  (H) 


(G-5a) 
(G-5b) 


which  indicates  that  equation  (G-4)  could  be  written  in  terms  of  the  - 
ion  with  equivalent  results,  as  required.  For  simplicity,  we  will  use 
equation  (G-5b)  to  write  A    in  terms  of  (H) ,   instead  of  (H) , , . 
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The  component  fluctuation  derivatives  in  terms  of  the  (H) . ,  are 


(compare  equations  4-55  to  4-57) 


N   ^%1 


kT  3N 


ol 


N   3%1 


(A),}   =  x^_(H)^/detU) 
—  11    o2 f—     — 


T,V,N 


o2 


kT  8N 


o2 


N    o2 


kT  9N 


T,V,N 


ol 


(A) 


-1 

12 


(A)"! 


ol 


T,V.N^2 


(G-6a) 


=  -X  _x  „(H),,/detU) 
ol  ol   —  1+     — 


(G-6b) 


N_!%2 

kT  9N  „ 
o2 


=  (A)  22 


[x  .+  x%(H)_J/detW) 
ol   ol  —  ii      — 


T,V,N 


ol 


(G-6c) 


where 


det(yl)  = 


x  ,x%(H),  +  x^x^J(H),  JH)^  -  (H)^ 


ol  o2  -'+-   ol  o2'-  -'11  -  +- 


1+-' 


(G-7) 


Double-Salt,  Single-Solvent  System 
Refering  to  equations  (4-72)  and  (G-3) ,  we  have 


A   = 


x  +  X  t(H).^ 
ol    ol  —  11 


X   X  -,(H) 
ol  o2  —  12+ 


X   X  ^(H) 
ol  o3  —  13+ 


X  X  „(H) 
ol  o2  —  12+ 


''ol''o3*--^13+ 


X 


— —  +  X  o(H)        X   X  t(H)   „ 
V      o2  —  2+2+    o2  o3  —  2+3+ 


^o2^o3^-^2+3+ 


+   X    t(H)^ , ^, 

V  o3  —  3+3+ 


(G-8) 


where,    from  equation    (G-la) , 


(H),2_   =    (H)i2+ 


(H)l3_   =    ©13+ 


(G-9a) 

(G-9b) 
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(H) 


+    (H) 


-'2+2-        V.   x_    ■    ^-'2-2-        v„.x7"^    ^-^2+2+ 


(H) 


2-  o2 

=        1 
3+3-        v_   X   „ 

3-  o3 


+    (H) 


2+  o2 

3-3-   =  ^^TTT""   ^^^3+3+ 
i+   oi 


(G-9c) 


(G-9d) 


©2+3+=    ^^>2-3+=    (li^2+3-   =    ^"^2-3- 


(G-9e) 


Rewriting^  in   terras   of    (H),,o    ,    and    (H) , , ^      give 


2+2-' 


3+3- 


A   = 


^ol"^  \1^^>11 


^01^02^^)12+ 


ol  o3  —  13+ 


^01^^02^12+ 


^2®2+2- 


"02^3(^)2+3+ 


"01^03(^)13+ 


"o2"o3(-)2+3+ 
"o3®3+3- 


(G-10) 


Finally,  the  component  fluctuation  properties  are  (compare  equations 
4-73  to  4-78) 


N 

^%1 

kT 

^^1 

2   2 


=  ^o2"o3f(^)  2+2-^^)3+3--  (H)2+3+]/detU)         (G-lla) 


^'^'"o2'^o3 


N   3%1 


kT  8N 


o2 


^ol^o2"o3f(-^)2+3+®13+-(^)l2+(^)3+3-^/^^tU)    (G-llb) 


T,V,N  ^,N  „ 
ol   o3 


N_!%2 

kT  8N  ^ 
o2 


2   2 


^'^'V'\3 


"01^03(^)3+3-  -^  "ol"o3f(^)  11(^)3+3-  -  (^)l3+J 


det(;l) 


(G-llc) 


N   3%2 


kT  SN 


o3 


T,V,N    N  - 
ol   o2 


"ol"o2"o3f(^)l2+(S)l3+  -  (li)2+3+(^)ll^ 

-  "01^02^03(^)2+3+  '    ^^^(^)  (G-lld) 


where 


de.M)  =  ^i>^^2<3((H)2^.,_(H)3,3_-(H)^^3J  +  ^l^Z^^os'®!!® 


2+2- 
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In  the  above  equations,  the  remaining  derivatives  of  the  chemical  poten- 
tial may  be  found  by  a  simple  exchange  of  subscripts  2  and  3. 

Double-Salt,  Single-Solvent  System  with  Common  Ion 
In  this  system,  A   can  be  determined  from  equations  (G-2)  and  (4-95) . 
Here,  A_   is  found  to  be  identical  to  that  of  the  previous  system  (equa- 
tion G-8) .   Now,  however,  the  z^   matrix  (equation  4-96)  differs,  leading 
to: 

©!_  =  (^'2-^o2+  ^3-\3)''f^2-^o2^«)l2++  ^3-^03^^^13+1  (G-12a) 

©24-  =  (^2-^o2+  ^3-^o3)"'[^  +  ^2-^o2(^>2+2++  ^3-^o3®2+3+l   (^'"^^b) 

(H)3+_  =  (V2_x^2+  ^3-^o3)"'f7^  ^  ^^2-^o2®  2+3++  V3_x^3  (11)3^3^]   (G-12c) 

2         2 

(H)    =  (v   x  +  V   x   )    [1  +  -^  X   +  -^  X   +  v^  X  9(H)  0.0^ 
I-   oz    3-  o3         v„    o2   \>  o3     2-  o2  —  2+2+ 

z+        J+ 

+  2v2_x^2^3_x^3(H)2^3^+  ^l.^l^^  ^^,^\  (G-12d) 

As  the  limit  to  one  salt  is  taken  (x  ^^0  or  x  _->-0)  ,  the  above  equations 

o2       o3 

agree  with  equations  (G-5) ,  as  expected.   However,  they  do  indicate  that 
the  interrelationship  between  the  (II )  .  ,  '  s  depends  on  the  system  of  in- 
terest ,  so  the  (]I)  .  .  cannot  be  uniquely  determined  for  two  ions  out- 
side of  the  framework  of  the  system.   Tliis  difficulty  again  makes  the 
direct  correlation  function  integrals  more  attractive,  despite  the  sim- 
plifications in  the  relationship  of  ^  to  H^  (which  are  still  much  more 
complex  than  that  of  <4  to  C  )  . 


APPENDIX  H 
INDEPENDENCE  OF  SALTS 


"Groups"  of  Salts 
The  simplest  method  to  assure  independence  of  the  salts  in  a  given 
electrolyte  solution  is  to  separate  the  salts  into  groups  having  the 
following  properties: 

1)  every  salt  in  a  given  group  has  an  ion  common  to  another 
salt  in  that  group, 

2)  no  salt  in  any  group  has  an  ion  common  to  any  salt  in  another 
group  (this  also  means  that  each  salt  is  in  only  one  group) , 

3)  the  groups  contain  as  few  salts  as  possible, 

4)  salts  with  no  ions  common  to  other  salts  form  single-salt 
groups. 

An  example  is  the  system  with  the  following  v  matrix  (ignoring  the 
solvents  which  form  their  own  groups),  with  seven  different  ions  (n=7) 

in  five  salts  (n  =5). 
o 


V  = 


V 


1,1+ 

0 
0 

^1,1- 
0 

0 

0 


V 


0 

2,2+ 
0 

0 

^2,2- 
0 

0 


0 

^3,2+ 
0 

0 

0 

^3,3- 
0 


V 


4,1+ 
0 

0 

0 

0 

0 


0 
0 

'5,3+ 

0 

0 
0 


^4,4-     ^,4-^ 
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We  group  the  \}_   matrix  according  to  the  above  rules  so  that  salts  1,4,5 
form  group  one  and  salts  2  and  3  form  group  two.   Salts  1  and  4  have 
a  common  cation,  4  and  5  a  common  anion  so  all  must  be  grouped  together. 
Salts  1  and  5  are  in  the  same  group  because  if  either  1,4  or  5  were 
placed  in  a  third  group,  rule  two  would  be  violated.   If  a  sixth  salt 
were  added,  having  no  ions  common  to  the  other  five  salts,  it  would 
form  a  third  group  by  rule  four.   The  most  complicated  systems  can 
easily  be  separated  into  groups  by  examining  their  v_  matrices. 

Tests  for  Independence 

After  the  salts  are  grouped  as  suggested  above,  each  group  can  be 
tested  for  independence  of  their  salts.   As  a  first  test,  if  the  number 
of  salts  in  a  group  is  greater  than  or  equal  to  the  number  of  ions  in 
that  group,  the  salts  are  not  independent.   For  systems  where  each  salt 
has  only  two  kinds  of  ions  (+  and  -) ,  each  group  must  have  exactly  one 
more  ion  than  salt  for  the  salts  to  be  independent. 

Allowing  for  salts  containing  more  than  two  kinds  of  ions  yields  a 
situation  that  is  somewhat  more  complicated.   In  this  case,  it  is  pos- 
sible to  have  fewer  salts  than  ions  and  still  have  dependent  salts.   A 
test  which  can  be  used  for  each  group  is  to  see  if  a  set  of  coefficients, 
E    ,  can  be  found  such  that,  in  the  sum  over  all  salts  in  the  group, 


y  /?  V.   =  0  for  all  i=l,...,n  (H-i) 

'^   a   la 
a 

Notice   that  H     =  0  if  salt  a  has  an  ion  that  is  not  common  to  any  other 

a 

salt.   Furthermore,  if  every  salt  has  such  a  unique  Ion,  the  salts  are 
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l»dei>endent,  but  this  is  not  a  necessary  condition.   For  example  for 
the  group  where 


V  = 


"1,1+ 
^1,2+ 

1^1,1- 


0 


2,2+ 


2,1-J 


we  see  that  salt  2  does  not  have  a  unique  ion,  but  salts  1  and  2  are 
independent. 


As  an  example  of  a  dependent  system  where  n=n  +1,  we  take  (n=5. 


n  =4) 
o 


V  = 


1,1+ 

0 

"3,1+ 

0 

1,2+ 

^2,2+ 

0 

0 

0 

^2,3+ 

0 

^,3+ 

1,1- 

^2,1- 

0 

0 

0 

0 

V.    ^ 

\>,  -. 

3,2-    4,2- 


Since  no  salt  has  a  unique  ion,  we  must  include  all  the  salts  in  the 
summation  of  (H-1) .   Using  the  charge  neutrality  conditions  on  each 
salt. 


y  V.  z.  =  0 


for  all  a=l , . . . ,n   , 


we  can  indeed  show  that  a  set  of  ff  's  exists  so  that  salts  are  dependent. 
In  general,  all  distinct  sets  of  R   's  must  be  found  (excluding  simple 
multiples)  and  for  all  such  sets,  one  salt  must  be  eliminated  as  de- 
pendent.  The  choice  for  the  salt  to  be  eliminated  can  be  found  among 
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thos'.i  N^'icb  have  no  unique  ions.   This  method  is  very  similar  to  that 
used  to  determine  independence  of  reactions  (Denbigh,  1971,  p.  169ff). 

After  each  group  of  salts  is  reduced  to  independent  salts  as  above, 
the  constraint  matrix,  z,  can  be  determined. 


APPENDIX  I 
ELECTROSTATIC  EXPRESSIONS  FOR  MULTISALT  SYSTEMS 


Debye-Hueckel  Limits 
The  formalism  of  Chapter  4  is  directly  applicable  to  a  general 
system  of  multiple  solvents  and  solvents.   The  Debye-Hueckel  terms  were 
given  only  for  single  solvent-single  salt  systems.   Since  we  deal  only 
with  electrostatic  effects,  these  terms  can  only  be  written  for  a  single 
solvent,  or  in  terms  of  a  "pseudo-component"  defined  as  the  salt-free 
mixture  of  the  solvents  and  designated  below  by  the  subscript  1.   For 
such  a  system,  the  activity  coefficient  limits  become 


n 

1        o        1 

V  2.nY.   =-Sq(p)xq)  /xix 
a   '±a      Y^a   ^  oy  r  (I-la) 

where 

2 

n 

q  =  y  V.  z.  (=0  for  solvent)  (I-lc) 

x=2 

These  summations  begin  at  2  since  we  include  only  the  salts,  not 
solvents.   The  partial  molar  volumes  become 

n 

_      _       T   1        °         1 

V  -  v"   =  ^  S  q  (p  y   X   q  )'^  (I-2a) 

oa    oa   2  v  a   ^„   oy  y 

Y=2 
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dine. 


(I-2b) 


Where  v   is  the  partial  molar  volume  of  salt  a  at  infinite  dilution  in 
oa 

the  solvent  mixture.   And  the  apparent  compressibility  becomes 


i-l<  '^Iw.''^ ^°^"..>' 


Y=2 


o  oy 
Y=2 


(I-3a) 


where 


S''"RT 


8ic, 


Y    r  2   .   „   "'^1 


s,  =  -V-{k:  +  2 


4   1     8P 


dine. 


-   6k 


1   8P 


+  9 


9«,ne, 


8P 


2    2 

-6  — ;r^ 


■'T    ap' 


}  (I-3b) 


The  Debye-Hueckel  limits  to  the  direct  correlation  function  integrals 
become: 


O  9 

1-C°  =  ^(1+  I  V  X  )  -  -^(l   V  X  ) 


-  {■ 


K  RT 


+ 


b     „   b,      ,      O 


3/2 


(I-4a) 


1-C 


V      „   S   3  £ne, 
og    3^   Y 1 


la   V  K,RT   2  V  K,   3P 
a  1        a  1 


q  p^  (  )  X   q  ) 
„^a  1    ^^  OY  Y 


Y=2 


(I-4b) 


and 


1-C 


1    h 


a6 


2v  v.(  y  X   q  )' 


(I-4c) 
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Mean  Spherical  Model  Expansion 
As  described  in  Chapter  5,  a  half-power  expansion  in  salt  density 
derived  from  the  Mean  Spherical  Model  is  used  to  improve  the  Debye- 
Hueckel  approximation  for  the  ion-ion  direct  correlation  function  inte- 
gral.  The  equations  for  these  terms  are  given  in  Chapter  5  also.   How- 
ever, actual  calculation  of  properties  requires  integration  of  these 
terms  over  composition  and  pressure.   The  composition  integrations  of 
the  Debye-Hueckei  terms  are  straightforward,  but  this  is  not  the  case 
for  the  expanded  terms  of  the  Mean  Spherical  Model.   Therefore,  the 
salt-density  expansion  of  the  activity  coefficient  of  salt  a,  including 
the  Debye-Hueckel  portion,  is  given  below 

V  any        =  -S^q  (p  I  x   q  )^+a^  f  v.  z^  ^  B   (o  n)K^       (1-5) 

2 
where  a    ,  n,  and  a     are  defined  in  Chapter  5,  and  k  is  the  inverse 
1   k 

Debye  length.   The  B   's  are  as  follows: 

^il  =  ^l^^i  +  \^ 

hi  =  ^2(2°i  +  H^   +|b2(a2+|n2) 

B.^  =  ha^ia^T)^)  +  ^^(o^r)^  +  a^n^  +  n^n2)  +  ^c^Co^  +  2n2) 

^14  "  ^4(\)(^°i  -  '-^^\^  +  b^(n^)(2a^n2  +  o^n^  +  i2n^n2) 


The  a  's,  b  's,  c.'s,  d,  and  e,  are  fitted  parameters  and  are  identical 
i     1     1     4      4 

to  those  used  for  the  direct  correlation  function  integrals. 


APPENDIX  J 
HARD-SPHERE  FORMULAE 


Hard-Sphere  Terms  for  Direct  Correlation  Function  Integrals 
The  hard-sphere  contribution  to  the  direct  correlation  function 
integrals  can  be  derived  from  the  mixture  form  of  the  Carnahan-Starling 
equation  of  state  (Mansoori,  et  al.,  1970).   For  two  species  (not 
components)  i  and  j  in  solution,  we  have 


,^o,.hs    (      .      )^ 


+  <;  (a.a.)^}/(l-u)^  +  9(0  a^y/a-t.   )^ 

CoCa.o  ) 
+  ^  ^  \  [9c-,(a.+o.)  +  6^,0.0,  +[6+^(-15+9^)]/c 
(1-Co)     ^  ^  J      ^  ^  J 


■^3 


where 


-  (o^+o  )C2[6+C3(-15+1253)]/C3  +  l^^o ^o ^\^^T,^(.-ir^r,^   x 


(26-14c3))]/C3(l-C3)}  +  (yi,  ^{o  ^o  ^'^  lx^{\-r,  ^ 


x{c3-(a.+a.)^2  +  C20io7c3}/C3  (-J-^^^ 


^."t  I' A  ''-'"' 
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Xlft©  c^'.tn't'iiciiL    functions   are,    of    course 


Unfortunately,  the  infinite  dilution  limits  of  these  correlation 
functions  are  not  correct,  so  the  model  must  be  altered  to  account  for 
this  deficiency.   This  accounting  is  given  in  Chapter  5,  where  the 
following  notation  is  used.   These  correction  terms  are  used  only  in  the 
solvent-solvent  and  solvent-ion  direct  correlation  function  integrals. 
If  we  allow  s  to  denote  the  solvent  pseudo-component  and  y  the  salts,  we 
have 


„hso     lim     ^hs 

^ss  =  (x   }.0  ^ss  (-^-23) 

oy 


,  ,     1  •  salts  ^  ^hs 

hs'  ^   lim  Y  d  C 

ss    {x   }^0  ^  dx 

oy  Y     oy 


(J-2b) 


As  a  work  of  explanation,  these  are  limits  as  the  total  salt  concentra- 
tion goes  to  zero  while  the  salt-free  solvent  composition  is  held 

constant.   We  also  define  C    similarly  to  C  ^°  in  equation  (J-2a) . 

sa  ss 

Finally,  the  hard-sphere  contribution  to  the  activity  coefficient  is 
given  by, 


iny      =   -£n(l-^)[l-3(-^)^  +  2i-^-^y] 

"■  T^  f Vi+3^2^i+3^1^'+3C3(^)'-C3(2-^3)(^)3] 
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1  r-^r    r   ^3^9   2^2      3   2,^2'^i,2, 

1  3   3  3    3 

+  3  [3daJ-r    d^f]  (J-3) 

(1-^3)  ^   "■        "^   ^   ^ 


Again,    this   is   for  a   species  and  not  a  component.      For  components,   we 
use 

n 

ilny     =      y   V.    toy.  (J-4) 

oa      .    -,    la        1 
1=1 

hs 
This  expression  is  the  analytic  integration  of  the  C   's,  which,  of 

course, for  our  calculations  are  integrated  from  a  reference  state 

r   r 
(T,  p  ,  X  ).   When  we  integrate  over  all  the  member  densities  to  T, 

p ,x  we  get  simply, 

Jiny   (T,p,x  )  -  Any    (T,p^,x'^) 
oa     -o      'oa      -o 


which  is  much  easier  to  evaluate  than  the  multiple  integrations  of  the 
C^^'s. 
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